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A new method for declaring erasures in a concatenated coding scheme is 
investigated. This method is used with the rate 1/2 K = 7 convolutional code and the 
(255, 223) Reed Solomon code. Errors and erasures Reed Solomon decoding is used. 
The erasure method proposed use a soft output Viterbi algorithm and information 
provided by decoded Reed Solomon codewords in a deinterleaving frame. The results 
show that a gain of 0.3 dB is possible using a minimum amount of decoding trials. 
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Chapter 1 


Introduction 


Concatenated coding systems are often used for forward error correction to obtain 
large coding gains when transmitting information over unreliable channels. One of the 
most popular concatenated coding systems is illustrated in Figure 1.1. This concatenated 
system uses a convolutional code with Viterbi decoding for the inner code, and a Reed 
Solomon code as the outer code. This is effective for a number of reasons. 
Convolutional codes provide sufficient random error correction, but tend to generate burst 
errors for low signal to noise ratio at the decoder output. Reed Solomon (RS) codes have 
significant burst error correcting capacity, but do not handle random errors very well. In 
this concatenated system, the inner convolutional code is used to correct the random 
errors, and although the Viterbi decoder will produce short burst errors at its output, the 
outer Reed Solomon code will be able to correct these bursts. The effects of these burst 
errors can further be reduced by using an interleaver between the inner and outer 
decoders. In addition, the Viterbi decoder can further improve the performance by 
accepting soft decisions from the receiver. 

The use of erasures is one way to increase the performance of Reed Solomon 
codes. Erasure decoding can be thought of as the simplest form of soft decision. An 
erasure indicates the reception of a signal whose corresponding symbol value is in doubt. 
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In some cases it is better to erase the symbol than to force a decision that may be 
incorrect. Erasing a position gives information to the decoder as to the location of a 
possible error. A block code with minimum distance d min can correct v errors and p 
erasures as long as the inequality 2 • v + p < d min is satisfied. Therefore, it is possible 

for a t-error correcting RS code to correct more than t errors if errors are transformed into 
erasures. 



Figure 1.1 A concatenated coding system using inner convolutional code with Viterbi 
decoding and an outer Reed Solomon code 

In the concatenated system in Figure 1.1, the Viterbi decoder produces hard 
outputs for input to the Reed Solomon decoder. The full capability of the concatenated 
system is not fully realized because no reliability information is exchanged between the 
inner and outer decoders. If the Viterbi decoder could be modified to generate reliability 
information about its output, this information could be used to declare erasures at the 
input to the Reed Solomon decoder, thus improving the performance. One method that 


2 












can be used to accomplish this is the Soft Output Viterbi Algorithm (SOVA). The 
method proposed by Hagenauer and Hoeher [4] uses information provided by the path 
metrics in the Viterbi decoder to determine a reliability value associated with each 
outgoing bit. 

One application where this gain could be potentially useful is in NASA deep 
space missions. The transmission of data over large distances, combined with limited 
transmission power, results in low signal to noise ratio at the receiving end. This, 
coupled with the fact that the data being transmitted is in the form of compressed images 
where the required probability of error is 10~ 5 , leads to the need for a powerful coding 
system [16]. The NASA standard for deep space communications is the (255, 223) 16 
error correcting RS code as the outer code, and the rate 1/2 convolutional code with 
constraint length K = 7. Interleaver depths of I = 2 to 8 have been used. The use of a 
SOVA and an errors and erasures RS decoder can provide additional gains with no need 
to modify the transmitting end. This enables erasure decoding to be used in existing 
missions. This is particularly helpful for missions where unforeseen problems occur. 
The Galelaio mission where the main antenna failed is one such instance. Every tenth of 
a decibel gain that can be obtained in this instance is extremely helpful [15]. 

One method used to improve the NASA standard for deep space communications 
through the use of erasures has been investigated by Paaske [7]. This method uses the 
deinterleaver to provide information concerning the probable locations of errors in non- 
decoded Reed Solomon codewords in an interleaving frame. In an deinterleaving frame 
there are I Reed Solomon codewords, where I is the interleaving depth. If after 
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attempting to decode the frame, some of the RS codewords fail to decode, redecoding is 
used. Erasures are declared using information provided by the error positions in the 
successfully decoded RS codewords. Because the Viterbi decoder produces burst of 
errors at it’s output, and the data is fed into the deinterleaver by row and output by 
column to the Reed Solomon decoder, the bursts occur at the same symbols in 
neighboring Reed Solomon words in the deinterleaver frame. If some but not all of the 
Reed Solomon words in the deinterleaving frame have been successfully decoded, the 
positions of the errors in the decoded words are known. The knowledge of the error 
positions can be used to declare erasures in the same positions in neighboring, yet to be 
decoded Reed Solomon codewords. 

1.1 Proposed research 

The purpose of this report is to investigate the performance of the use of a Soft 
Output Viterbi Algorithm used in a concatenated coding scheme with an errors and 
erasures RS decoder. The reliability information provided by the SOVA will be 
converted into Reed Solomon symbol erasures for the RS decoder. A table of least 
reliable symbols will be compiled for each RS codeword, and systematically erased. In 
addition, another method loosely based upon Paaske’s method will be investigated. This 
method combines the SOVA output with a deinterleaver. The table of least reliable 
symbols can be modified using additional information provided by the deinterleaver. If 
after the first decoding of a deinterleaving frame, there are less than I successful decoding 
of RS codewords, redecoding is attempted. It turns out that not only does the SOVA 
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output produce burst errors, but the reliabilities for these error symbols are identical. This 
information is used to modify the table of smallest reliabilities. The performance of these 
codes will be obtained through the use of a computer simulation written in C computer 
language. The convolutional code developed for use in this simulation is capable of 
handling any code rate and constraint length. The Reed Solomon code, likewise, can 
handle any symbol size and number of symbol errors corrected. The Reed Solomon 
decoder is an errors and erasures decoder. Although the codes developed are capable of 
handling any size code, the NASA standard coding system will be investigated with 
various interleaving depths. The simulation will be performed over a AWGN channel 
using BPSK modulation and Raised Cosine FIR filters. The coding systems will also be 
simulated over an ideal BPSK channel. 

The structure of this report is as follows. Chapter 2 contains all of the background 
information. Chapter 3 will contain the details of the computer simulation. Chapter 4 
will contain the strategy for declaring RS symbol erasures from the reliability information 
generated by the SOVA, and the strategy for using the SOVA with the deinterleaver for 
redecoding. Chapter 4 will also present the results of the simulation for both methods 
investigated. Chapter 5 will contain conclusions, and ideas for possible future research. 
The simulation flow charts are found in Appendix A and the C language source code used 
to perform the simulations can be found in Appendix B. 
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Chapter 2 


Background 

Before discussing the two methods for erasure declaration presented in this report, 
it is helpful to become familiar with some of the basic concepts of error control codes. 
This chapter will contain all of the background necessary to understand the various 
elements used in the concatenated system. Encoding and decoding of Reed Solomon and 
convolutional codes will be reviewed. In addition, the method used for errors and 
erasures decoding in the Reed Solomon decoder will be discussed, as well as the method 
used for generating the soft outputs in the Viterbi decoder. Block interleaving will be 
briefly discussed, in addition to the redecoding method proposed by Paaske. 

2.1 Reed Solomon codes 

Bose-Chadhuri-Hocquenghem (BCH) codes are a powerful class of cyclic codes 
which outperform all other block codes with the same block length and code length [9]. 
These codes are a generalization of Hamming codes to allow multiple error correction. 
Reed Solomon (RS) codes are special subclass of BCH codes which utilize non-binary 
symbols. The non-binary symbols used in RS codes are formed using finite field 
arithmetic. Finite fields are sometimes called Galois fields and are denoted by GF(p), 
where p is the number of elements in the field, and is a prime number. 
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2.1.1 Galois fields 


A field is a set of elements in which we can do addition, subtraction, 
multiplication, and division without leaving the set. Subtraction and division are defined 
by the additive inverse and the multiplicative inverse. Addition and multiplication must 
also satisfy the commutative, associative, and distributive laws. A field with a finite 
number of elements is called a finite field. For example, GF(7) = {0, 1, 2, 3, 4, 5, 6} is a 
field under modulo 7 addition and multiplication. Reed Solomon codes are codes with 

symbols from the field GF(2 m ) , where GF(2 m ) = {0, l,a,a 2 ,...,a 2m ~ 2 }. The field 

GF(2 m ) is an extension of the ground field GF(2), and the elements in this field can be 
represented by an ordered sequence of m components, (a 0 , a,, a 2 , ..., a m l ), or an m- 

tuple. Each of the components are from the ground field GF(2). The 2 m elements of the 
field GF(2 m ) are defined by an irreducible polynomial, or a primitive polynomial P(x). 
Each element will satisfy the condition P(a) = 0. The primitive polynomials that define 
the elements for m = 3 to 9 are as in Table 2.1. For example, the elements in GF(2 3 ) are 
defined using the primitive polynomial 1 + x + x 3 . The elements either are the zero 
element ‘O’, the identity element ‘1’, or some power of the base element a. The element 
a 3 is derived from the primitive polynomial and the relationship P(a) = 0. 

P(a) = a 3 + a + 1 = 0 
or a 3 = a + 1 . 

All other elements are simply generated by multiplication by a. The table repeats after 
a m 2 (i.e. a 6 • a = 1 , a 6 a 2 = a etc.). The elements for GF(2 3 ) are as follows. 
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0 


0 


1 = 

1 


a = 

la = 

a 

a 2 = 

a • a = 

a 2 

a 3 = 

1 + a 


a 4 = 

a • a 3 = 

a(a + 1)= a + a 2 

a 5 = 

4 

a • a = 

a(a + a 2 ) = 1 + a + a 

a 6 = 

a • a 5 = 

a-(l + a + a 2 ) = 1 + a 


Table 2. 1 : List of primitive polynomials for m = 3 to 9 


m 

P(X) 

3 

1 +x + x 3 

4 

1 +x + x 4 

5 

1 + X 2 + X 5 

6 

1 + X + X 6 

7 

.1 + x 3 + X 7 

8 

1 + x 2 + x 3 + x 4 + x 8 

9 

1 + x 4 +x 9 


It is useful to represent these elements in a number of ways. The polynomial 
representation is given by a 0 + a,a + a 2 a 2 + ... + a m . 1 a m_1 and the m-tuple 




representation is given by (a 0 , a 15 a 2 , a m .,). The elements and the various 

representations for GF(8) are given in Table 2.2. 


Table 2.2 Three representations for the elements of GF(8) 
generated by 1 + x + x 3 


Power 

Polynomial 

3-tuple 

Representation 

Representation 

Representation 

0 

0 

(0 0 0) 

1 

1 

(10 0) 

a 

a 

(0 10) 

a 2 

a 2 

(0 0 1) 

a 3 

1 + a 

(1 10) 

a 4 

a + a 2 

(0 11) 

a 5 

1 + a + a 2 

(1 1 1) 

a 6 

1 + a 2 

(10 1) 


Multiplication and addition follow the rules of finite field algebra. Multiplication 
of two elements is accomplished by adding the powers of the two elements modulo 2 m -1. 
For example, in GF(8), a 4 - a 5 = a (4+5)mod7 = a 2 . For addition of two elements in a 
field, it is useful to use the m-tuple representation of an element. Consider a = 
(a 0 , a,, a 2 , ..., a m .,) and b = (b () , b , , b 2 , ... , b m .,) . The addition of a and b is simply 

the addition of each component in the m-tuple representation, 
namely (a 0 + b 0 ,a, + b, , a 2 + b 2 , ..., a ra _, + b m .,). Because each component of the 
m-tuple is from GF(2), binary addition is used. 
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2.1.2 Generating Reed Solomon codes 


A t-error correcting Reed Solomon code with symbols from GF(2 m ) has the 
following parameters: 

Block Length n = 2 m - 1 

Number of information symbols k = n - 2t 
Minimum Distance d min = 2t + 1 

The generator polynomial of a t-error correcting Reed Solomon code is: 
g(x) =(x + a)-(x + a 2 )---(x + a 2t ) 

where g(x) has all of its roots and coefficients from GF(2 m ). The code generated from 
g(x) is a (n, n - 2t) cyclic code. The code words are generated by: 

c = QG 

Where G is the generator matrix in systematic form. Let us design a t = 2 error correcting 
Reed Solomon code using symbols from GF(2 3 ) = GF(8). We know that: 

Block length n = 2 3 - 1 n = 7 

Information symbols k = n - 2t k = 3 

The generator polynomial for this (7, 3) Reed Solomon code is given by: 

g(x) = (x + a) • (x + a 2 ) • (x + a 3 ) • (x + a 4 ) 
or 

g(x) = a 3 + a • x + x 2 + a 3 • x 3 + x 4 
= [ a 3 a 1 a 3 1] 
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The generator matrix in non-systematic form is: 


G(X) = 


a' 

0 a 


a 1 a 

3 


1 0 0 


a 


1 a" 1 0 


0 0 a 3 


a 


1 a 1 


To get the matrix into systematic form, we must convert the last three columns into an 
identity matrix. This is accomplished by row operations. The result after doing so is: 



V 

a 

1 

a 3 

1 

0 

0 

G(X) = 

a 6 

a 6 

1 

a 2 

0 

1 

0 


a 5 

a 4 

1 

a 4 

0 

0 

1 


The information bits to be transmitted are u = [ 010 Oil 110 ]. From Table 2.2 we 
know that these bits correspond to the symbols u = [a a 5 a 3 ] in GF(8). Using 
c = u-G, we obtain the code vector 




a 3 

a 

1 

a 3 

1 

0 

o’ 

c= [a 

a 3 a 3 ]• 

a 6 

a 6 

1 

a 2 

0 

1 

0 



a 5 

a 4 

1 

a 4 

0 

0 

1 

c = 

[ a a’ 

a 4 

a 4 

a 

a 5 

a 

3 ] 



c = [010 110 011 011 010 011 110] 


The encoding of Reed Solomon codes can be also accomplished using a shift 
register circuit. For a t-error correcting RS code, the generator polynomial is given by: 
g(x) =(x + a)-(x + a 2 )---(x + a 2t ) 


n 



= g» + g|X + g,X 2 + ... + g 2l .,x 2, -‘ + X 2 ' 

where g(x) has all of its roots and has coefficients from GF(2 m ). The generator 
polynomial g(x) has been chosen so that it and codewords generated by it have zeros for 
2 • t consecutive powers of a 

g(a') = 0 for j = 1, 2, ... 2 t 

The code generated from g(x) is a (n, n - 2t) cyclic code. The encoding of a non-binary 
cyclic code is similar to the encoding of a binary cyclic code. Let 

u(x) = u 0 + u j x + u 2 x 2 + ••• + u 2M x 2t_l 

be the message to be encoded. In systematic form, the 2t parity check symbols are the 
coefficients of the remainder b(x) = b 0 + b,x + b 2 x 2 + ••• + b 2M x 2t_1 which is 

obtained by dividing the message polynomial u(x) by the generator polynomial g(x). In 
hardware, this is accomplished by using the shift register circuit of Figure 2.1. The 
encoder circuit works as follows. The k information symbols are first loaded into the 
circuit. At the same time, the k information symbols are transferred directly to the output. 



Figure 2. 1 Encoding circuit for t error correcting RS code 
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When all of the information symbols have been read in, the 2t parity symbols are present 
in the 2t registers denoted b 0 , b,, b 2M , and are then transferred to the output, thus 

completing the systematic code word. This process can best be illustrated with an 
example. Consider the t = 2 error correcting (7, 3) Reed Solomon code. The generator 
polynomial is 

g(x) = (x + a) • (x + a 2 ) • (x + a 3 ) • (x + a 4 ) 

= a 3 + a ■ x + x 2 + a 3 • x 3 + x 4 

and it’s corresponding encoding circuit is given in Figure 2.2. The encoding of the 
information symbols u = [a a 3 a 3 ] is given in Table 2.3. The information symbols 



Parity Check Digits 

Figure 2.2 Encoding circuit for t = 2 error correcting (7, 3) RS code 

are fed directly to the output to the encoding circuit. At t = 1, the first information 
symbol is fed into the encoding circuit, and the register contents are modified. 
Information symbols are fed into the encoder until t = 3. At this time, the 4 parity 
symbols are present in the registers b 0 , b p b 2 , b 3 , and are sent to the output of the 
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encoding circuit. The encoded vector is equal to c = [b 0 b 1 b 2 b 3 u 2 u 1 u 0 ] or 
c = [ a a 3 a 4 a 4 a a 5 a 3 ]. 


Table 2.3 Shift register contents for the encoding of u = [a a 5 a 3 ] 


t 

Input Symbols 

Gate 

b o 

b, 

b 2 

b 3 

0 

- 

- 

0 

0 

0 

0 

1 

c 

o 

II 

a 3 

a 6 

a 4 

a 3 

a 6 

2 

u, = a 5 

a 5 + a 6 = a 

a 4 

1 

a 2 

a 6 

3 

u 2 = a 

a + a 6 = a 5 

a 

a 3 

a 4 

a 4 


2.1.3 Decoding of Reed Solomon Codes 

Let r(x) = r 0 + r,x + ... + r n .,x n ~‘ be a received polynomial which is equal to 
a codeword c(x) = c 0 + c,x + ... + c^x" -1 corrupted by an error pattern 
e(x) = e 0 + e,x + ... + e n .,x n_1 . 

r(x) = c(x) + e(x) 

The syndrome of the received polynomial is obtained by evaluating r(x) at the 2 • t zeros. 
Sj = r(a J ) = c(a j ) + e(a j ) j = 1, 2, ..., 2 • t 

Any codeword c(x) will have zeros for these 2 • t powers of a, and thus, have a syndrome 

equal to zero. Therefore, the syndrome of the received word is equivalent to the syndrome 
of the error pattern. 
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( 2 . 1 ) 


S j = r(a J ) = e(a j ) = ]T e k *(a j ) k j= 1,2,..., 2 t 


k=0 


If there are v errors in positions i , , i 2 , . . . , i v ,(2.1) can be expressed as 


S, = Xe, ,(«')' 


1 = 1 


= e ; a 1 + e ; a 2 + ... + e ; a 


(2.2) 


j = 1, 2, ..., 2 • t 


To reduce the notational complexity of (2.2), the error locations will be defined as 
X, = a M , and the error magnitudes as Y, = e i] where 1=1,2, ..., v. (2.2) then becomes 


Sj = YjXj + Y 2 X 2 + ... + Y V X V 

S 2 = YjXj 2 + Y 2 X 2 2 + ... + Y V X V 2 

s 3 = ^Xj 3 + Y 2 X 2 3 + ... + Y V X V 3 (2.3) 


S 2t = YjXj 21 + Y 2 X 2 2t + ... + Y v X v 2t 
The error locator polynomial A(x) is defined as 


A(x) = (1 - xX,)(l - xX 2 ) (1 - xX v ) 

= 1 + A,x + A 2 x 2 + ... + A v . 1 x v 1 


+ Ax' 


(2.4) 


where the roots of A(x) are the error locations X,, X 2 , ..., X v . The coefficients of the 
error location polynomial A, 1 = 0, 1, ..., v are related to the error locations by the 
following equations 


A 0 = 1 


+ X., , + X. 


A, — Xj + X 9 + 

A 2 = X,X 2 + X,X, + 

A 3 = X,X 2 Xj + X,X 2 X 4 + ... + X V . 3 X # ,X V + X V . 2 X V ,X, 


+ X V . 2 X V + x..,x 


V-l V 


(2.5) 


A v = X,X 2 X 3 ••• X V .,X, 
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(2.3) and (2.5) are related by Newton’s identities [11] 

Sj + Aj = 0 
S2 + AjSj + 2A2 = 0 


S 3 + A,S 2 + A 2 SJ 

+ 3A 3 = 0 

S v + AjSy.j + A 2 S v _ 2 + ... 

+ A v .|Sj + 

S v+i + AiS v + A 2 S v _! + 

... + A v S t 

S 2t + AjS 2t _i + A 2 S 2t _ 2 + 

... + A v S 2t . 


( 2 . 6 ) 


(2.6) can be solved directly to obtain the coefficients of the error locator polynomial, but 
such methods require a number of computations proportional to t 3 [2]. This makes a 
direct solution of (2.6) not practical, especially for RS codes that need to correct a large 
amount of errors. Berlekamp’s algorithm is much more computationally efficient method 
of correcting RS codes. The complexity increases linearly with t, so codes correcting 
large numbers of errors can be implemented [11]. Berlekamp’s algorithm first finds a 
minimum degree polynomial A (1) (x) whose coefficients satisfy Newton’s first identity. 
This polynomial is tested whether the second Newton identity is also satisfied. If it does, 
then A ( 2 ) (x) = A (l) (x) . If not, then a correction term, or discrepancy is added to A (1) (x) 
to form A ( 2 ) (x) such that A (2) (x) satisfies the first two Newton’s identities. Next 
A (2) (x) is tested whether it satisfies the third Newton’s identity, etc. This process 
continues until A (2t) (x) is obtained. Then A(x) = A (2t) (x) . If there are less than t errors, 
A(x) produces the error pattern. 
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Massey’s shift register based interpretation of Berlekamp’s algorithm is known as 
the Berlekamp-Massey algorithm [2, 11]. The Newton’s Identities in (2.6) can be 
expressed in an alternate form 

V 

Sj = - X A i s j-i j = V + I - v + 2 ’ •••> 2v - (2.7) 

i = l 

Massey [17] recognized that (2.7) can be represented physically using a linear feedback 
shift register (LFSR) as shown in Figure 2.3. 



Figure 2.3 LFSR interpretation of (2.7) 


The output of the LFSR will be the 2t syndromes S 1? S 2 , S 2t , and the register taps 
are the coefficients of the error correction polynomial A(x). The LFSR can be designed 
to generate the known sequence of syndromes such that A(x) is of the smallest degree. 
The procedure for finding the taps of the LFSR is similar to Berlekamp’s algorithm. First 
a connection polynomial T(x) = 1 + AjX + ... + A l .jX L_1 + A L x L is formed 
whose coefficients are the taps of a length L LFSR. The Berlekamp-Massey algorithm 
first finds T(x) of length L = 1 such that the first output of the LFSR is the first syndrome 
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Sj . The second output of the LFSR is compared to the second syndrome, and if the two 
are not equal then the connection polynomial is modified using a discrepancy term. If the 
two are equal the taps remain the same. The third output of the LFSR is compared to the 
third syndrome, and if they are not equal, the taps of the connection polynomial are 
modified. This process continues for 2t iterations. At the end of the 2t iterations, the taps 
of the LFSR specify the coefficients of the error correction polynomial A(x). The details 
of the algorithm are presented below. 

2.1.3.1 Berlekamp-Massey Algorithm [2] 

1. Compute the syndrome of the received codeword S } = r(a J ) . 

2. Initialize the following variables 

Error locator polynomial A(x) = 1 4 

Index r = 0 

Temporary storage B(x) = 1 
Shift register length L = 0 

3. Set r = r + 1. 

4. Compute the r th discrepancy, which is the error in the next syndrome 

A r = i A j S r-j 

j=0 

5. If A r = 0 , set B(x) = x • B(x) and go to step 11. 

6. Compute the new connection polynomial 

T(x) = A(x) - A r • x • B(x) 
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7. If 2 L > r - 1 , set B(x) = x • B(x) and go to step 10. 

8. Store old shift register after normalizing 

B(x) = A r _l A(x) 

9. Update shift register length L = r - L. 

10. Update the shift register 

A(x) = T(x) 

11. If r < 2 • t , go to step 3. 

12. If deg A(x) * L , there are more than t errors. Stop. 

13. Determine the roots of A(x). The inverses of these roots are the error locations 

X„x 2 X v . 

14. Determine the corresponding error values Y,, Y 2 , ..., Y v . 

The simplest method to find the roots of A(x) in step 13 is by using a process 

known as a Chien Search. This is a trial and error approach which computes A(a j ) for j 

= 0, 1, .. 2 m -2. If A(a j ) = 0 , then a j is a root of A(x). The error magnitudes can be 

calculated by using the Forney algorithm[3, 11]. First, compute a syndrome polynomial 
S(x) from the 2 • t syndromes. 

S(x) = 1 + (2.7) 

j= 1 

The error evaluator polynomial Q(x) can be computed by the product of the syndrome 
polynomial S(x) and the error locator polynomial A(x). 

Q(x) = S(x)-A(x) modx 2t+l (2.8) 
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Next, compute the derivative of the error locator polynomial A'(x) . 


A'OO - Ix.nd - xXj) 


(2.9) 


i = 1 j * i 

The error magnitudes Y, , Y 2 , . . . , Y v can then be calculated using the Forney Algorithm. 


Y, 


x.-'ncx,- 1 ) 

A'cxr 1 ) 


1=1,2, v 


(2.10) 


For example, let r(x) = a + a 2 x + a 4 x 2 + a 4 x 3 + a 6 x 4 + a 5 x 5 + a 3 x 6 be a 

code word corrupted by an error pattern e(x). The first step in decoding is to compute the 
syndrome of the received polynomial. 

Sj = r(a) = a + a 2 a + a 4 a 2 + a 4 a 3 + a 6 a 4 + a 5 a 5 + a 3 a 6 
S, = 1 

S 2 = r(a 2 ) = a + a 2 a 2 + a 4 a 4 + a 4 a 6 + a 6 a 8 + a 5 a 10 + a 3 a 12 

5 2 = a 2 

5 3 = r(a 3 ) = a + a 2 a 3 + a 4 a 6 + a 4 a 9 + a 6 a 12 + a 5 a 15 + a 3 a 18 

5 3 = 1 

5 4 = r(a 4 ) = a + a 2 a 4 + a 4 a 8 + a 4 a 12 + a 6 a 16 + a 5 a 20 + a 3 a 24 
S 4 = a 6 


The syndrome polynomial is 

2-t 

S(x) = 1 + X s j xi = 1 + x + a 2 x 2 + x 3 + a 6 x 4 . (2.11) 

j = i 

Next, the Berlekamp-Massey algorithm is used to find the error locator polynomial A(x). 
The results of the computations for each iteration of the algorithm are given in Table 2.4. 
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Table 2.4 Results of the computations for each iteration of the 
Berlekamp-Massey algorithm 


r 

A, 

T(x) 

B(x) 

A(x) 

L 

0 

- 

- 

1 

1 

0 

1 

1 

1 + X 

1 

1 + X 

1 

2 

a 6 

1 + a 2 x 

X 

1 + a 2 x 

1 

3 

a 5 

1 + a 2 x + a 5 x 2 

2 4 

a + ax 

1 + a _ x + ax 

2 

4 

0 

1 + a 2 x + a 5 x 2 

a 2 + a 4 x 

1 + a 2 x + a 5 x 2 

2 


The error locator polynomial is found to be A(x) = 1 + a z x + a 5 x 2 . The roots of 


A(x) are a 6 and a 3 , and the inverses of these roots give the error locations X. = a and 


X 2 = a 4 . The error evaluator polynomial is 


£2(x) = S(x)A(x) = (1 + x + a 2 x 2 + x 3 + a 6 x 4 )(l + a 2 x + a 5 x 2 ) 


5 „ 2 


= 1 + ax + ax 


and the derivative of the error location polynomial is A'(x) = a 2 . The error values can 


be calculated using (2.10) 

a(l + a 6 a 6 + a 5 a 12 ) 


Y, = 


Y 2 = 


a' 


a(l + a + a ) 5 

2 = a 

a 


a 4 (l + a 6 a 3 + a 5 a 6 ) a 4 (l + a 2 + a 4 ) 


= a 


a‘ 


a‘ 
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The error polynomial is e(x) = a 5 x + a 5 x 4 , and the corrected polynomial is 

r(x) + e(x) = (a + a 2 x + a 4 x 2 + a 4 x 3 + a 6 x 4 + a 5 x 5 + aV) + (a 5 x + a 5 x 4 ) 

= a + a 3 x + a 4 x 2 + a 4 x 3 + ax 4 + a 5 x 5 + a 3 x 6 

2.1.3.2 Errors and erasures RS decoder 

In order to correct both errors and erasures, certain modifications to the 
Berlekamp-Massey algorithm will need to be made. Suppose that a received polynomial 
r(x) contains v errors in locations i p i 2 , i v and p erasures in locations 

j p j 2 , j p . An errors and erasures Reed Solomon code can correct v errors and p 
erasures as long as 2-p + v < d min where d min is the minimum distance of the code. 
The error locations are given by X k = a' k k = 1, 2, ..., v and the erasure locations are 
given by U, = a Jl 1 = 1, 2, ..., p. The erased positions are known at the beginning of 
the decoding operation, and are filled with zeros before the decoding begins. 

2.1.3.3 Berlekamp-Massey algorithm for errors and erasures [2] 

1 . Substitute zeros into the erased positions in the received word. 

2. Compute the syndrome of the received codeword Sj = r(a j ) . 

3. Initialize the following variables 

Errors and erasures locator polynomial A(x) = 1 
Index r = 0 

Temporary Storage B(x) = 1 


22 



Shift register length L = 0 

4. Set r = r + 1. 

5. If r > p go to step 10 

6. A(x) = A(x)-(1 - U r • x) 

7. B(x) = A(x) 

8. L = L + 1 

9. Go to step 4 

10. Compute the r th discrepancy, which is the error in the next syndrome 

L 

A r = Z A j S r-j 

j = 0 

1 1. If A r = 0 , set B(x) = x • B(x) and go to step 11. 

12. Compute the new connection polynomial T(x) = A(x) - A r • x • B(x) 

13. If 2 • L > r + p - 1 , set B(x) = x • B(x) and go to step 16. 

14. Store old shift register after normalizing 

B(x) = A r -'A(x) 

15. Update shift register length L = r-L + p. 

16. Update the shift register 

A(x) = T(x) 

17. If r < 2 • t , go to step 4. 

18. If deg A(x) ^ L, 2p + v > d mn . Stop. 
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19. Determine the roots of A(x). The inverses of these roots give the error locations 
X p X 2 , X v and the erasure locations U, , U 2 , U . 

20. Compute the error evaluator polynomial 


£2(x) = S(x) • A(x) mod x 


2t+l 


21. Use the Forney Algorithm to compute the error values 


X k fl(X k -') 

A'(X k ->) 


Y; 


and the erasure values 


i = 1, 2, ..., v 


(2.13) 


u.nqjr 1 ) . . , n 

Jl a'cu,- 1 ) J ’ ’■"’ p 


(2.14) 


For example, let r(x) = a + a 3 x + f x 2 + a 4 x 3 + a 6 x 4 + fx 5 + a 3 x 6 be a code 
word corrupted by an unknown error pattern e(x) and a erasure pattern 
f(x) = f x 2 + f x 5 with known positions and unknown values denoted by T . The first 

step in decoding is to insert zeros into the erased positions and compute the syndrome of 
the received polynomial. 


S, = r(a) = a + a 3 a + a 4 a 3 + a 6 a 4 + a 3 a 6 
S, = a 

S 2 = r(a 2 ) = a + a 3 a 2 + a 4 a 6 + a 6 a 8 + a 3 a 12 
S 2 = a 


s 3 = r(a 3 ) = a + a 3 a 3 + a 4 a 9 + a 6 a 12 + a 3 a 18 

5 3 = a 6 

5 4 = r(a 4 ) = a + a 3 a 4 + a 4 a 12 + a 6 a 16 + a 3 a 24 
S 4 = 0 


24 



The syndrome polynomial is 


s( x ) = i + X^j x ' = 1 + ax + a 6 x 2 + a 6 * 3 ■ 

i = i 

Next, the modified Berlekamp-Massey algorithm is used to find the error locator 
polynomial A(x). The contents of the variables for each iteration of the algorithm are 
given in Table 2 . 5 . 


Table 2.5 Results of the computations for each iteration of the errors 
and erasures Berlekamp-Massey algorithm 


r 

A r 

T(x) 

B(x) 

A(x) 

L 

0 

- 

- 

1 

1 

0 

1 

- 

- 

1 + a 2 x 

1 + a 2 x 

1 

2 

- 

- 

1 + a 3 x + x 2 

1 + a 3 x + x 2 

2 

3 

a 3 

1 + a 2 x 2 + a 3 x 3 

4 4 2 

a + x + ax 

1 + a 2 x 2 + a 3 x 3 

3 

4 

a 2 

1 + a 6 x + a 4 x 3 

4 2 4 

a x + x + a x 

1 + a 6 x + a 4 x 3 

3 


The error locator polynomial is found to be A(x) = 1 + a 6 x + a 4 x 3 . The roots of 
A(x) are a 2 , a 3 , and a 5 . The inverses of these roots give the error location X, = a 4 

and the erasure locations U, = a 2 and U 2 = a 5 , which were known at the beginning 
of the decoding operation. The error connection polynomial is 

Q(x) = S(x)A(x) = (1 + ax + a 6 x 2 + a 6 x 3 )(1 + a 6 x + a 4 x 3 ) 
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and the derivative of the error location polynomial is A'(x) = a 6 + a 4 x 2 . The error 
value can be calculated using (2. 13) 


v _ 0^(1 + a 5 oc 3 + ora b + orx y ) a 4 (l + a + a + a 4 ) 

Yj - 


, 2_6 


, 2 „ 9 


a 6 + a 4 a 6 


a 6 + a 3 


= a 


and the erasure values by using (2.14) 


Y, = 


or(l + q 5 q 5 + ora 10 + a 2 x 15 ) a 2 (l + a 3 + a 5 + a 3 ) 


a 6 + a 4 a i0 


a + 1 


= a 


Y - a + a 5 ** 2 + + a 2 x 6 ) 


2 4 


2 „ 6 


a 6 + a 4 a 4 


a (1 + 1 + a + a) « 

= a- 

a + a 


The error and erasure polynomials are e(x) = a 5 x 3 and f(x) = a 4 x 2 + a 5 x 5 and the 
corrected polynomial is 

r(x) + e(x) + f(x) = (a + a 3 x + a 4 x 3 + a 6 x 4 + a 3 x 6 ) + (a 5 x 3 ) + (a 4 x 2 + a 5 x 5 ) 

= a + a 3 x + a 4 x 2 + a 4 x 3 + a 6 x 4 + a 5 x 5 + a 3 x 6 

2.2 Block Interleaving 

Interleaving is commonly used to break up correlated errors into random errors by 
rearranging the symbols. This is done because most block and convolutional codes are 
optimal for random errors. Interleaving causes correlated errors to be spread out over 
time, and then the coding system can handle the errors as if they were random. There are 
two major types of interleaving, block and convolutional. Block interleavers are used in 
conjunction with the concatenated systems in Figure 1.1, and will be the only method 



discussed here. Because the interleaver is to be used in conjunction with the symbol 
based Reed Solomon encoder and decoder, the interleaving will be done on a symbol 
level, rather than on a bit level. Block interleavers can be implemented using an MxN 
matrix. The symbols are fed into the matrix by column, and fed out by rows. At the 
deinterleaving stage, the symbols are fed into the matrix by row, and output by column. 
Consider this simple example. The sequence {0, 1, 2, 3, ..., 11} will be fed into a 3x4 
block interleaver by column 

0 3 6 9 

1 4 7 10 

2 5 8 11 

The interleaver then outputs the data by row. The output sequence is {0, 3, 6, 9, 1, 4, 7, 
10, 2, 5, 8, 11}. The deinterleaving is accomplished by entering the sequence by row, and 
outputting by column. The deinterleaved sequence is {0, 1,2, 3, ..., 11}. 

2.2.1 Redecoding of deinterleaving frame using erasure 

Paaske [7] has developed a strategy to declare Reed Solomon symbol erasures 
using information provided by the deinterleaver. In order to understand the erasure 
declaring procedure discussed later in the report, and for comparison purposes, an 
overview of Paaske’ s method is presented below. 

Because the Viterbi decoder produces burst errors at it’s output, and the 
deinterleaver spreads these bursts over several codewords, it is highly likely that the burst 
errors will occur at identical positions in neighboring Reed Solomon codewords (RSW) 
in each deinterleaving frame. The output of the Viterbi decoder is fed into the 
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deinterleaver by row, and then the columns are fed to the RS decoder as in Figure 2.4. 
Each column of the deinterleaver makes up a Reed Solomon codeword. Assume after 
errors only decoding of each of the RSW, that some of the codewords in the deinterleaver 
frame were decoded correctly ( < 16 errors in the codeword) and some were undecodable 
( > 16 errors in the codeword). Because a RSW with v errors and p erasures can be 
corrected if 2«v + p < d min , declaring erasures and redecoding the deinterleaving 

frame may provide improvement. This improvement is highly dependent upon if the 
declared erasure hits an error. Erasures that hit errors will be called good erasures (GE) 
and ones that do not hit an error will be called bad erasures (BE). 


From Viterbi decoder 



Figure 2.4 Typical deinterleaving frame 
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Paaske’s method for declaring erasures in the de-interleaving frame makes use of 
the bursty nature of the Viterbi decoder output. Because the data is fed into the 
deinterleaver by row, the bursts will span over many RSW as in Figure 2.5. Let RSW(i) 
denote the i th codeword in an interleaving frame. It should be noted that a burst of 
length 1 at symbol k starting at RSW( ij ) will affect symbol k in RSW( i j +j) for ij + j < I 
and symbol k + 1 in RSW( i } +j-I) for ij + j > I . If codeword RSW(i) has been correctly 
decoded, the positions where errors have occurred will be known. It is highly probable 
that the same symbols in the neighboring codewords will also be errors. Paaske 
developed 4 erasure declaring procedures (EP1-EP4), three of which declare erasures in 
the non-decoded RSW using information provided by the decoded RSW in the 
deinterleaving frame. A brief description of the procedures is presented below. 

Symbol 



= Symbol error 


Figure 2.5 Typical burst error in a deinterleaving frame 

EP1: Assume that two RSW have been decoded, and both contain errors in 

position k. Let RSW( ij ) and RSW( i 2 ) be two correctly decoded RSW with an error in 
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position k in both as shown in Figure 2.6. Also assume that RSW(ij +1) to RSW(i 2 -l) 
have not been decoded. It is highly probable that position k is an error in these 
undecoded words, and is erased. These erasures are classified as double sided erasures 
(DSE) and the probability that a DSE is a GE is 0.96 [7]. An example of a DSE is shown 
in Figure 2.6. 

EP2: Assume that a codeword RSW(i) has been successfully decoded and that 
RSW(i - 1) and RSW(i + 1) have not been decoded. For all error positions in RSW(i) 
erase the same positions in RSW(i - 1) and RSW(i + 1). These are called single sided 
erasures (SSE) and are GE with a probability of 0.60 [7]. An example of SSE declaration 
is given in Figure 2.7 where RSW(i) is the decoded word with errors in positions k t , k 2 , 
k 3 , and k 4 . Note that the SSE declared in RSW(i + 1) at position k 2 and RSW(i - 1) at 
position k 3 do not hit symbol errors, and are therefore BE. 



RSW 

= Symbol Error DSE = Double sided E = Error corrected ir 

erasure decoded word 


Figure 2.6 An example of a double sided erasure declaration 
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EP3: Assume that RSW(i) has been decoded and contains e errors, and RSW(i - 1) 
and RSW(i + 1) have not been decoded. Also assume that S[ DSE can be obtained if EP1 
is used. The s { DSE are combined with s 2 SSE chosen from the e - s t possible SSE. 
The optimal choice for the number of SSE s 2 is treated in [7]. 


Symbol 

Number 



RSW 


= Symbol Error SSE = Single sided E = Error corrected in 
erasure decoded word 


Figure 2.7 An example of single sided erasure declaration 
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EP4: This procedure assumes that one of the non-decoded RSW has 17 errors. 
Two symbols are selected and make erasures. If there are 17 errors, the probability that 
an erasure is a good erasure is 1/15. 

The erasure declaring procedure proposed by Paaske [7] involves the following 

steps: 

1) Try decoding each RSW in the deinterleaving frame using errors only decoding. 

2) Set i d equal to the number of successfully decoded RSW. 

3) If i d = I, go to step 9. 

4) If i d = 0, go to step 8. 

5) Attempt to decode each non-decoded RSW using EP1 . 

6) Attempt to decode each non-decoded RSW using EP2. 

7) Attempt to decode each non-decoded RSW using EP3. 

8) Attempt to decode each non-decoded RSW using EP4. 

9) Stop. 

For steps 5) through 8), decoding is attempted on the first non-decoded RSW 
using the EP specified in the step. If this decoding attempt is successful, then proceed to 
step 2). If not successful, then try the next non-decoded RSW using the EP specified in 
the step. This continues until either one of the non-decoded RSW is successfully 
decoded, or all non-decoded RSW have been tried and none are successful. If all non- 
decoded RSW have been attempted using the given EP, and there are no successfully 
decoded RSW, then proceed to the next step. Erasure procedures EP3 and EP4 involve 
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selecting erasures in a systematic way. This step is repeated on each codeword until 
either a successful decoding of the codeword, or a maximum number of trials T max has 
been attempted. In the simulations conducted by Paaske, T max = 500 trials. 

2.3 Convolutional codes 

Convolutional codes are fundamentally different than block codes. Block codes 
divide the information sequence into segments of length k, and map these k bits onto a 
codeword of length n. Convolutional codes on the other hand convert the entire data 
stream into one code word, regardless of the length of the information sequence. A (n, k, 
m) convolutional encoder has k inputs and n outputs, where k < n and both k and n are 
small integers. The memory order m should be made large in order to achieve a high 
degree of error correcting capability [6]. 

2.3.1 Convolutional encoder 

Convolutional codes are implemented using a linear feed forward shift register 
circuit. A typical encoding circuit is given in Figure 2.8, and will be used as a model for 
discussing convolutional encoders and the Viterbi decoder. The information sequence 
u = (u 0 , u, y u 2 , ... ) is fed into encoder circuit k bits at a time. The memory elements 
are tapped and the bits contained in memory are added together using modulo-2 adders to 
obtain a pair of output data streams v (0) = (v 0 <0) , Vj (0) , v 2 (0) , ...) and 

v (l) = (v 0 (1) , Vj (1) , v 2 (1) , ...). These output sequences are combined to create the final 
codeword v = (v 0 <0) , v 0 (,) , v/%/ 0 , v 2 (0) , v 2 (1) , ...). 
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The constraint length K is defined as the maximum number of output bits that can 
be affected by any input bit. Since each information bit stays in the encoder for m + 1 
time units, and during each time interval the encoder produces n output bits, the 
constraint length is defined as K = n • (m + 1) . 



Figure 2.8 A (2, 1, 2) convolutional encoding circuit 


The structure of convolutional encoders can be expressed in a number of ways. 
One of these ways is using the impulse response of the encoder. The impulse response of 
the encoding circuit is obtained by letting the input u = (1 0 0 0...) and observing the 
output sequences as u enters the encoding circuit. An encoder with memory m generates 
an impulse response of length m + 1 . The impulse response, also known as the generator 
sequence, is written in the form g (0) = (g 0 (0) , g, (0) , g 2 (0) , ... , g m (0> ) and 

g (1) = (g 0 (1) > gi 0) > g 2 0) ’ ••• » g m (I) )- F °r the encoder in Figure 2.8 the generator 
sequences are g (0) = (10 1) and g (!) = (1 1 1). The two encoder output sequences can 
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be thought of as a linear convolution of the information sequence with the impulse 
response. The encoding equations can be written as 

v<») = u * g (0) (2.15a) 

v"’ =u*g l,) (2.15b) 

where * denotes discrete convolution using modulo-2 operations. The output at time t = t 
can be written as 

m 

V T U) = Z U T-i ' gi j 

i=0 

= u t g 0 (J> + u t _,-g, (J) + u,_ 2 -g 2 (JI + ... + u,_ m - g m a) . (2.16) 
For the encoder of Figure 2.5, (2.15) reduces to 

v t ,0) = u x © u t . a (2.17a) 

v,' 01 = u t © u t ., © u,. 2 (2.17b) 

where 0 denotes modulo-2 addition. The encoding of the information sequence 
u = (101101) is illustrated in Table 2.6. At time = 0, the contents of memory are 
initially set to zero. At time = 1 , the first information bit is fed into the encoder, and the 
output of the encoder is obtained by using (2.17). This process continues until all of the 
information bits have entered into the encoder. At this point, information bits are still 
contained in memory. Two more clock cycles are needed to move the last bits through 
the encoder. ( k • m ) zeros are fed into the input to move the last information bits through 
the encoding circuit. The encoded sequence is v = (1 1 01 00 10 10 00 01 11). 
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2.3.2 Decoding of Convolutional codes 

The Viterbi algorithm is a widely used method for the decoding of convolutional 
codes. The algorithm was developed by A. J. Viterbi in 1967 [6], and is a maximum 
likelihood decoder. 


Table 2.6 Encoding of the information sequence u = ( 1 0 1 1 0 1 ) 


Time 

Input 

m, 

m 2 

v (0) 

v (,) 

0 

- 

0 

0 

- 

- 

1 

1 

1 

0 

1 

1 

2 

0 

0 

1 

0 

1 

3 

1 

1 

0 

0 

0 

4 

1 

1 

1 

1 

0 

5 

0 

0 

1 

1 

0 

6 

1 

1 

0 

0 

0 

7 

- 

0 

1 

0 

1 

8 

- 

0 

0 

1 

1 


2.3.2.1 State Diagram 

For every encoding circuit, there will be a corresponding state diagram. An 
encoder with memory m will have 2 m possible states, illustrating the contents of the shift 
registers in the encoding circuit. There are k binary inputs to the encoder for each clock 

cycle, which results in 2 k branches entering and exiting each state. The process details 
used to create the state diagram are given in Table 2.7. For the encoding circuit of Figure 
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2.8, there will be 4 states, and 2 branches entering and exiting each state. The state 
diagram for the circuit in Figure 2.8 is given in Figure 2.9. 


Table 2.7 Development of the state diagram 


Initial State 

m, 

m 2 

Input u 

Output v 

New State 

S () 

0 

0 

0 

00 

So 

So 

0 

0 

1 

1 1 

s, 

S, 

1 

0 

0 

0 1 

s 2 

S, 

1 

0 

1 

1 0 

S 3 

S 2 

0 

m 

0 

1 1 

So 

S 2 

0 

i 

1 

00 

S, 

S 3 

1 

i 

0 

1 0 

s 2 

s 3 

1 

i 

1 

0 1 

S 3 


0 / 00 



Figure 2.9 State transition diagram for the encoder in Figure 2.8 
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2.3.2.2 Trellis Diagram 

A trellis diagram is a state diagram extended to include the passage of time. The 
encoding of a sequence of data corresponds to a unique path through the trellis diagram. 
The trellis diagram for the state diagram in Figure 2.8 is shown in Figure 2.10. If an 
information sequence of length k • L bits is fed through an encoding circuit, where L is 
the total number of k bit codewords, the resulting codeword will be of length 
N = n • (L + m) bits. Each of the 2 L code words of length N is represented by a unique 
path through the trellis. 



Figure 2.10 Trellis diagram (based on the encoder in Figure 2.8) 


The convolutional coding problem is shown in Figure 2.11. Assume an information 
sequence u = (u 0 , Uj, ... , u L _j) of length k-L bits is encoded into a code word 
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y = (y 0 . y I » ••• * y L +m-i) of length N = n-(L + m) bits. A noise-corrupted version of 


the transmitted sequence r is received where r and y have the following form 


_ ( r (0) (!) (n-1) (0) (1) (n-1 ) \ 

' A 0 > A 0 > ••• ’ x 0 » M ’ *1 ’ ••• ’ 1 L+m-l ) 



Noise 


Figure 2.1 1 Convolutional coding system 


The Viterbi algorithm generates the estimate y' of the transmitted sequence r which 
maximizes the conditional probability p(r I y) . Assuming the channel is memoryless, 

each received bit will be independent of the noise process affecting all of the other bits. 
Therefore p(r I y) can be expressed as 

L+m- l 

p(rly) = fJtpO^I y i (0) )-p(r i (1) ly i (1) ) ••• p(r j (n ' 1) ly. (n " ,) )] 

i=0 

L+m- l n-j 

= 11(11 p( r i <j) 1 yi lj) ) (2.i8) 

i=0 j=0 

The log-likelihood function is obtained by taking the logarithm of each side of (2. 18). 

L+m-l n-I 

logp(rly) = Y, X lo g P( r i U) I yi (j) ) (2.19) 

i=0 j=0 
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This is done because it is, in general, easier to implement summations rather than 
multiplication in hardware. The log likelihood function, logp(rly), is called the Path 

Metric associated with the path y and is denoted M(r I y) . The terms log pCr^ ly i (j) ) are 
called Bit Metrics 

M(rj 0) ly^) = log p(r i ci) I y i (i> ) (2.20) 

In the hardware implementation of the Viterbi decoder, it is more convenient to use 
positive integers for the metric values rather than the actual bit metrics. This can be 
accomplished by using 

M(r i ®ly j ^ ) )= a-[ log p(r i < J ) ly,^') + b) (2.21) 

where a and b are chosen to obtain small, positive integer values for the metrics which 
can be implemented easier in hardware. The path metric for a codeword y is then 
calculated as 

L+m-I n-1 

M ( r| y)= X X M ( r ®'y ®)- (2.22) 

i=0 j=0 

The k th branch metric for a codeword y is defined as the sum of the bit metrics 

M (f k I y k ) = XM(r k (j) I y k (j) ). (2.23) 

j=o 

The k th partial metric for a path is obtained by summing all of the branch metrics for the 
first k branches the path follows. 

k-i 

M k (rl y) = X M( >i ly,) (2.24) 

i =o 

= X £ M(r®ly i*)- (2.25) 

i=0 j=0 
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The Viterbi Algorithm finds the path through the trellis with the largest path metric, 
which is the maximum likelihood estimate y' of the received word r. At each time 
interval the algorithm, computes the partial metrics entering each state. The largest 
metric is chosen as the surviving path at each state, and all other paths entering that state 
are discarded. This process is continued until the end of the trellis is reached. The final 
surviving path is the maximum likelihood estimate y of the codeword. 

2.3.2.3 The Viterbi Algorithm [6] 

1. At time t = m, compute the partial metric for the single path entering each 
state. Store the value of this metric at each state. 

2. Increase t by 1. Compute the partial metric for the path entering each state. 
This will be equal to the branch metric entering the state plus the surviving metric from 
the previous state. Out of the 2 k paths entering each state, the path with the largest metric 
is chosen and the remaining paths are discarded. The metric of the surviving path is 
stored at each state. 

3. If t < L + m, repeat step 2. If not, stop. At time t = L + m, all paths have 
returned to the all zero state. There will be only one path remaining, and this path is the 
maximum likelihood estimate y' . 

2.3.2.4 Hard Decision Decoding 

In hard decision decoding, the receiver determines whether a zero or one was 
transmitted. These zeros and ones are the input to the Viterbi decoder. If the channel is 
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memoryless and if the probability of a bit error is independent of the transmitted bit, then 
the channel is said to be a binary symmetric channel (BSC). The BSC is shown in Figure 
2.12 where p is the probability that a bit is in error. The conditional probabilities for the 
BSC are given in Table 2.8. 



1 

Re ceived 
Symbol 

0 


Figure 2. 12 Binary symmetric channel model 


Table 2.8 Conditional probabilities for BSC 


p ( r i <j) i yi U> ) 

o 

II 

II 

>> 

j-s 

II 

O 

1-p 

p 

T® = 1 

p 

1-p 


For the BSC, choosing a = [log 2 p- log 2 (l - p)]“‘ and b = -Iog 2 (l-p) in 
(2.21) yields the bit metrics in Table 2.9 [6]. The maximizing of the bit metrics 
M(r i (i) I y, 0 ) coincides with the minimization of the Hamming distance. For the BSC 
case, the path metric is simply the Hamming distance d(r, y), and the Viterbi algorithm 
will choose the surviving paths as the ones that have the minimum partial path metrics. 
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Table 2.9 Metric table for BSC 


M( !•« 1 y,' 1 ’) 

o 

n 

y/ j > = 1 

r = 0 

0 

1 

r = 1 

1 

0 


Consider the information sequence u = (101 101) that was encoded using the (2, 1, 
2) convolutional encoder in Figure 2.8. The encoded sequence is y = (1 1 01 00 10 10 00 
01 11). If this sequence is transmitted on a BSC, and no errors occur in the transmission, 
r = y. The decoding of this received sequence is illustrated in Figure 2.13. Note that the 
final path has a path metric value of 0. The decoded sequence is ( 1 0 1 10 1) which is 
obtained by tracing back the maximum likelihood path noting the input bit associated 
with each branch. 


2 3 0 2 3 



Figure 2.13 Hard decision Viterbi decoding 
of r = (11 01 00 10 10 00 01 11) 
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If the same sequence y = (1 1 01 00 10 10 00 01 1 1) is transmitted and the received 
sequence is r = (11 1 1 00 10 00 00 1 1 1 1) where the erroneous bits are denoted with 
the bar over the bit, the decoding is illustrated in Figure 2. 14. The decoded sequence is ( 1 
0 1 10 1), which is identical to the information sequence. The decoder corrected the 
three errors in the received sequence. 



Figure 2.14 Hard decision Viterbi decoding 
of r= (11 I TOO 10 00 00 Tl 11) 

2.3.2.5 Soft Decision Decoding 

Hard decision decoding simply assigns a zero or a one at the receiver, utilizing 
only two decision regions. Soft decision decoding makes use of q-bit quantization which 
results in multiple decision regions ranging from a “strong-one” to a “strong-zero”. 
Using soft decisions results in approximately 2 dB gain over hard decision Viterbi 
decoding [9]. A discrete memoryless channel (DMC) is shown in Figure 2.15. The 
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DMC is completely described by a set of transition probabilities between a zero or a one 


being transmitted and one of Q = 2 q levels at the receiving end. It has been found that 



Received 

Symbol 


Figure 2. 15 DMC channel model for Q = 4 levels 


using 8-level quantization results in only a 0.25 dB loss when compared with using 
infinitely fine quantization [11]. Consider the transition probabilities given in Table 2.10. 
The modified metric table is given in Table 2.1 1, and is obtained by using (2.21) with b 
= 1 and a = 17.3. In choosing a and b, a is typically chosen to obtain a metric value equal 
to zero for the smallest metric value. The metric values obtained by using (2.20) are, in 
general, real valued. Simply rounding these values off to the nearest integer may lead to 
round off errors. The scaling factor b in (2.21) is chosen to make the metrics as close as 
possible to being integer values, while keeping the values as low as possible. This will 
reduce some of the error that may occur when rounding off the metric values. 
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If the same codeword y = (11 01 00 10 10 00 01 11) is transmitted over the 
channel, and the received sequence is r= (1,1 2 0,1, 0,0 2 1 2 1 2 1 2 0, 0,0 2 0 2 0 2 0 2 0,), 

then the decoding process using soft decision Viterbi decoding is illustrated in Figure 
2.16. The decoded sequence is u = ( 1 0 1 10 1), which is identical to the information 
sequence that was transmitted. 


Table 2.10 Conditional probabilities for DMC 


p ( r i (j) 1 yj (J) ) 

o 

J 1 

y ; (j) = 1 

r i <j) = 0, 

.4 

.1 

r i W = 0 2 

.3 

.2 

= 1 2 

.2 

.3 

r i <j> = 1, 

.1 

.4 


Table 2. 1 1 Metric Table for DMC 


M( r, (j) 1 y, <J) ) 

O 

II 

y t (i> = i 

r = 0, 

10 

0 

r = 0 2 

8 

5 

r = 1 2 

5 

8 

r = 1, 

0 

10 


2.3.2.6 Truncation length 

In practice, information sequences are very long. It is not practical to wait until 
the entire sequence is received to begin decoding. This would result in long delays and 
require large amounts of storage. It has been found that a decision can be made on the k 
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information bits that were received (t - 5) time units before, where 5 is called the decision 


depth, or truncation length. If the truncation length is made 4 to 5 times the constraint 
length, there will be very little loss in performance [10]. 


18 33 69 74 92 



0,1, 0 t 0 2 1 2 1 2 1,0, 0,0, o 2 o 2 0,1, 


Figure 2.16 Soft decision Viterbi decoding of 
r = (1,1 2 0,1, 0,0 2 1 2 1 2 1 2 0, 0,0 2 0 2 0 2 0 2 0,) 

The implementation of the truncated Viterbi decoder makes use of 2 m shift 
registers, each of length k -8. At any time t, there are 2 m surviving paths, with one 
surviving path terminating in each of the 2 m states. For each surviving path, the only 
information that needs to be stored are the information bits associated with that path. No 
information about the route the path took is necessary, just the information (output) bits 
associated with that path. At time t, n bits are input into the decoder. The branch metric 
is calculated, and the surviving path is chosen as the path with the largest metric, as in the 
standard Viterbi decoder. The path information for each state at time t is equal to the 
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previous state at time (t - 1) shifted k bits to the left. The k information bits associated 
with the surviving branch at time t are then shifted into the register. Consider the Viterbi 
soft decoding example given in Figure 2.16. For a decoding depth of 5 = 5, the shift 
registers, surviving branches, and metric values are shown in Figure 2.17. At time t = 0, 
the shift registers are empty. At each time interval, the survivor branch is chosen, and the 
contents of the shift register at the pervious state are copied to the new state, shifted to the 
left, and the information bit associated with the branch is inserted into the register. This 
process continues until the register is full (i. e. t > 8). At this point, the path with the 
highest metric is chosen as the surviving path. Only one bit is output at a time, so this 
corresponds to the leftmost position in the shift register, which was the information bit 
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t = 0 t = 1 t = 2 1 = 3 t = 4 t = 5 t = 6 


Figure 2. 17 Shift register contents for the soft decision decoding 

in Figure 2.16 
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inserted t - 5 time units before. In Figure 2.17 at t = 5 = 6, the path with the highest 
metric has a path metric of 87, and terminates in state S 2 . The leftmost bit in the shift 
register is equal to 1, and is the output for this time interval. Now, since t > 6, the 
decoder can take in n bits, compute the path metric and determine the surviving path for 
each state. The decoder then chooses the path with the maximum metric, and outputs the 
leftmost bit. At time t = 6 in Figure 2.17, the path with the highest metric has a path 
metric equal to 105, and terminates in state Sj . The leftmost bit in that shift register is 0, 
and is the output bit. This process continues for the remainder of the decoding operation. 


2.3.3 Soft Output Viterbi Algorithm 

The Viterbi algorithm can be modified to give either a reliability value or a 
probability that a given bit is correct. The method used to implement the Soft Output 
Viterbi Algorithm (SOVA) is based upon Hagenauer and Hoeher’s method [4], For 
simplicity, this discussion will only consider convolutional codes where k = 1. The 
reliability of a binary random variable can be defined in terms of a log-likelihood value 


L(u) = log 


Prob(u = 1) 
Pr ob(u = 0) 


(2.26) 


where the sign of L(u) corresponds to the hard decision (i.e. if L(u) > 0, u = 1 and if L(u) 
< 0, u = 0) and the magnitude IL(u)l is the reliability of this decision. The larger the 
magnitude, the greater the reliability of the decision. 
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At each of the 2 m states at time t, the Viterbi decoder selects the survivor path as 
the one with the largest path metric. The accumulated path metric and the path 
probability are related by 

t n-l 

M(r I y) = log 2 p(rly) = J I logp(ij (j) I yj (j) )• (2.27) 

i=l j=o 


If the bit metrics are given by M(r i (j) ly i ( J ) ) = a-[log 2 p(r i ( J ) ly^) + b], then the 
path probability is 


p(r I y) = 


M(r I y) - ntb 

e a 


(2.28) 


Each state S k (k = 0, 1, ..., 2 m -l) will have two entering paths, a survivor path with 
metric M, and a competing path with metric M 2 (M, > M 2 ). The probability of 
choosing the wrong path is given by 

Prob(path 2) 

pe = — - (2 29) 

k Prob(path 2) + Prob(path 1) 


——M2 - nkb 
e a 


— — 'M 2 - nkb - nkb 

a + e a 


(2.30) 


1 In 2 

= t- where A = (M, - M,) (2.31) 

1 + e a 

With this probability the Viterbi decoder has made an error in the path it has chosen as 
the survivor path. Consider the two paths merging in state S 0 at time t in Figure 2.18. 

The all zero path is the survivor path with metric M, , and the other is the competing path 
with metric M 2 . The two paths are the same up to time t - 5 m ( 5 m = 6 in this case). At 
this point, the paths diverge and there are three differing information bits between the two 
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0 



M2 


Figure 2.18 Example of the SOV A 

paths at times t-2, t-3, and t-5. Using the probability of selecting the wrong path given in 
(2.31), the probability associated with each bit pj can be modified for all bits Uj for 

times j = t, t-1, ..., t - 8 m [4]. 

Pj = Pj if u,. = u 2j (2.32a) 

Pj = Pj(l-Ps t ) + (l-Pj)Ps k if u ij 54 u 2j (2.32b) 

where u } . and u 2 . are the output bits at time j on paths 1 and 2 respectively. The first 

case can be neglected because choosing path 2 instead of path 1 would result in no error 
for the j th bit. Because the case in (3.32a) can be neglected, there is only a need to check 
u ij ^ u 2 j f° r ti me s t-m, t-m-1, ..., t-8 m +l. This is because in order to terminate in any 

given state at time t, the m input bits prior to time t must be equal in order to create the 
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given state, S k . The probability in (2.32b) can be transformed into a bit reliability. From 
(2.26), the log-likelihood for the j th bit can be expressed as 

L J = log (2.33) 

Pj 

This can be combined with (2.31) and (2.32) to obtain an expression for updating the 
likelihood function. 


Lj = log 


1 - PjO-Psk) + U - Pj )Psk 
PjO-Psk) + (l-Pj)Psk 


= log 


1 • pj( TT^ ) ' ( 1 ' p i ) rr^ 


= log 


Pj + e A - Pje A 

Pje A + 1 - Pj 


1 " Pj 
1 + 11 


= log- 


Pi 


A 1 ' Pj 
e + - 


P 


Lj = log 


< L; + A 
1 + e J 

A L : 
e + e J 


(2.34) 


A good approximation of this expression is to simply take the minimum of L } and A as 


the new reliability [4]. 


Lj = min(Lj, A) (2.35) 

For register exchange mode with truncation depth 6, each state S k will have a shift 
register of 8-q bits where q-1 bits represent the magnitude of the likelihood value Lj 
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and one bit for the sign of Lj , which corresponds to the output bit Uj . The procedure 

can be summarized as follows. For every state at time t, compute the path metrics for the 
two paths entering state S k . Choose the path with the higher metric as the survivor path 

and update the path information in the register for S k . The reliability at time t, L t is 

In ^ 

initially set to «>, Compute the metric difference A = -(M, - M 2 ) . For j = t - m, t - 

a 

m - 1, t - 5 m compare the information of the two paths. If Uj * u 2 , then update the 

new reliability using (2.35). After the surviving paths in each state have been determined, 
and the reliability information has been updated for each state, the state with the highest 
path metric is determined and the reliability and output bit at t - 8 are the output for the 
decoder for time t. 
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Chapter 3 


Simulation Techniques 

Computer simulation plays an important role in the design and testing of 
communications systems [5]. The results obtained from simulation can give a good 
indication of how an actual system will perform under a variety of conditions. 
Performance evaluation of complex communication systems using analytical techniques 
can be difficult, if not impossible. Testing of the concatenated coding schemes presented 
in this report using computer simulation will give an estimate of how these codes perform 
under realistic conditions. In addition, the codes will be simulated using an “ideal” 
AWGN channel to test how the coding systems will perform under ideal conditions. 
Monte Carlo techniques were used to obtain the results for both simulations. Monte 
Carlo techniques are relatively simple to implement. Data is generated at the input to the 
simulation. This data is then run through the system being simulated. The data at the 
output of the simulation is compared to the data at the input to determine the number of 
errors. The probability of an error is simply the number of errors divided by the total 
number of simulation points. To be statistically confident in the results, the simulation 
should produce at least 50 errors [5]. For small values of bit error rate, a large number of 
bits will be needed, which results in longer processing times. This large amount of 
processing time is one of the drawbacks of using Monte Carlo techniques. The first step 
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in simulating a communications system is to describe the system in block diagram form, 
where each block represents a signal processing operation. The model used for this 
simulation is given in Figure 3.1. 



Simulation can either be done at baseband or bandpass. Baseband simulations 
have no carrier frequency. This reduces the complexity of the system models such as the 
filters. Bandpass simulations require a higher sampling frequency than baseband, and 
therefore, more computational time. Bandpass simulations are necessary when studying 
upconversions, downconversions, and the effects of adjacent channel interference. Most 
simulations that involve a single information signal can be done at baseband [13]. The 
simulation used in this report was performed at baseband. 

3.1 Random Number Generators 
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The various signals that exist in communication systems are random in nature. 
The information signals found in communication systems typically use random voltage or 
current waveforms to transmit information from one place to another. Noise is an 
unwanted random signal, and causes errors in the information being sent. In order to 
represent the random signals found in communication systems, a random number 
generator (RNG) will be needed. Random number generators do not produce truly 
random numbers, but a sequence of “pseudo” random numbers which repeat after a 
period of time. These sequences should be stationary and uncorrelated. The mean, 
variance, and other parameters computed for different segments of the RNG sequence 
should be equal for a RNG to be considered stationary [5]. Having a RNG with a period 
less than the simulation length will cause correlation in the RNG sequence. Choosing a 
RNG with a large period is desirable to avoid correlation. 

3.1.1 Uniform Random Number Generator 

Uniform RNG generate equiprobbile numbers within a given interval, typically 
between zero and one. Uniform RNG can be generated using the following multiplicative 
congruential algorithm 

. Ij+j = a • I j (modulus m) (3.1) 

where a and m are integer constants. If a and m are chosen carefully, (3.1) will produce a 
sequence of random numbers with a maximum period of m. The random number 
generator used in this simulation can be found in [8] and has a very long period of 

* 2.3 x 10 . This is accomplished by combining two RNG sequences with 
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m 1 = 2147483563 
aj = 40014 
m 2 = 2147483399 
a 2 = 40692. 

In addition to a long period, this RNG passes all of the relevant statistical tests [8]. 
Assuming the probability of a ‘O’ and a ‘1’ are equal, the random bit stream can be 
obtained by using a uniform random number with the following conversion 
bj = 0 if Ij < 0.5 

bj = 1 if Ij > 0.5 

3.1.2 Gaussian Random Number Generator 

There are situations that call for random numbers with different distributions, and 
are typically generated by performing a transform of a uniform deviate. The generation of 
Additive Gaussian White Noise (AWGN) in simulations calls for a sequence of normally 
distributed random numbers. This normally distributed sequence can be generated by 
using the Box Muller Method [5, 8]. If Xj and X 2 are two independent variables with 
uniform distribution between 0 and 1 , then 

Y 1 = jll + a • yj - 2 • ln(Xj ) • cos(2 • n • X 2 ) (3.2a) 

and 

Y 2 = jll + a • 2 • ln(Xj ) • sin(2 • n • X 2 ) (3.2b) 

are independent Gaussian variables with mean ji and standard deviation a. 
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3.2 Sampling 

In order to represent the signal in the simulation, the signals will need to be 
sampled. The Nyquist rate of 2 • B is the minimum sampling rate for bandlimited signals 
of bandwidth B. For simulations, the sampling rate needs to be much higher to accurately 
represent the analog signal and to reduce the affects of aliasing. The number of samples 
per symbol should be an even integer between 8 and 16. Having the number of samples 
per symbol greater than 16 is not necessary for most simulations [13]. 

3.3 Filters 

In Communication systems, filters are needed for the purposes of bandlimiting 
signals and rejecting out of band noise. The filtering in the process can produce 
something known as intersymbol interference. 

3.3.1 Intersymbol Interference 

Consider the effects of passing a series of impulses spaced T b seconds apart 

through a low pass filter. Each impulse produces its own output from the filter. The 
output from one pulse extends into the output of pulse that starts T b seconds later. This 

is known as Intersymbol Interference (ISI) and it can produce errors at the receiver. 

The effects of ISI can be avoided by an appropriate choice of a low pass filter. 
Nyquist proposed a technique that theoretically produces zero ISI. This is accomplished 
by creating in the receiver a pulse that resembles the sin x/x shape, crossing the axis at 
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intervals of T 5 . The receiver samples the incoming wave at intervals of T b , so at the 
sampling instant, the tails of the preceding outputs are crossing the axis, and are therefore 
zero. The only non-zero component is the pulse to be sampled, which solves the problem 
of interference from other symbols. The filter proposed by Nyquist is the “Raised 
Cosine” filter and theoretically produces zero 1ST The transfer function for this filter is 


H(f) = 


1.0 


for |f| < fj 


cos2 ( I7 r — 4\ ) for f > < l f l - 2 B o- f i 

4 • (rJ 0 - tj ) 

0 for |f| > 2 • B 0 - fj 


R 

where B 0 = is the filter bandwidth, a is the rolloff factor. The frequency fj and 


the bandwidth Bq are related by 

a = 1 - — 0 < a < 1 

B 0 

The frequency response for different rolloff factors is given in Figure 3.2. The minimum 

Rl 

bandwidth value of B n = — - is obtained when a = 0.0. This value of rolloff is not 

u 2 


obtainable in practice. Practical filters use rolloff values ranging from 0.2 to 1.0 [ 12 ]. 

In some applications, the filtering operation to produce zero ISI needs to be split 
between two filters, with one at the transmitter and one at the receiver. The optimum 
partition in the sense of optimizing the signal to noise ratio is to divide the filter transfer 
function equally between transmit and receive filter [5, 18]. 

H Transmit (0 = H Rgceive (f) = -y/H(f) 
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This is known as a Square Root Raised Cosine filter. 


H(f) 



f 


Figure 3.2: Frequency responses for different rolloff factors 


The Raised Cosine filter produces zero ISI only when driven by an impulse. If the 
filter is not driven by an impulse, then the transfer function of the filter must be divided 
by the Fourier transform of the input signal. For a NRZ (Non-Return Zero) square pulse 
train, the Fourier transform has a spectrum with a sin x/x shape. The transfer function 
then becomes 


; rc-f-T, 

sin(7t • f • T s ) 


for |f| < fj 


H(f) = 


n-f-T s 
sin(7i • f • T s ) 


• cos 2 


, n-(f-fi) 

MBo-fj) 


) 


0 


for fj < |f| < 2-B 0 - fj 
for |f| > 2 • B 0 -fj 


(3.3) 
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where T s is the symbol period. Note that the raised cosine filter is bandlimited to 


• (1 + a) , and at f = R s , x/sin x goes to infinity. Therefore, for this system to work, a 
must be less than 1.0 [12]. 

3.3.2 Digital Filters 

In this simulation, filters will be needed for the purpose of band limiting signals, 
producing zero intersymbol interference (ISI), and rejecting out of band noise. Because 
this simulation is in discrete time, digital filters will be needed to accomplish the above 
objectives. There are two major types of digital filter design methods: Infinite Impulse 
Response (HR), and Finite Impulse Response (FIR). Both have advantages and 
disadvantages, and neither is best for all situations. The optimal filter design method 
must be determined by analyzing the requirements and objectives of the application. 

The HR method uses the widely available filter functions from analog filter 
theory. This method starts with an analog filter transfer function, and then translates this 
analog function in such a way that makes it suitable for discrete-time systems. Filters 
designed using the HR method will be recursive in nature (the output of the filter depends 
on previous filter outputs, as well as past and current values of the input), and the filter’s 
impulse response will be infinite. HR filters require fewer coefficients than FIR filters, 
and have a closed form design technique that does not require iteration. Some of the 
disadvantages of HR filters include non-linear phase response, and the use of feedback in 
the implementation that can cause instability if not carefully implemented. 
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The FIR method does not rely on analog filter theory. Instead the frequency 
response of the desired filter is used to determine the digital filter coefficients. This 
design method is non-recursive in nature, and the impulse response has a finite number of 
terms. Filters designed using the FIR method are always stable, and have a linear phase. 
FIR filters need a high number of coefficients to adequately describe the impulse response 
of the filter. This large number of filter coefficients results in longer processing times, 
and can be a great disadvantage if used in real time applications. In addition, the design 
procedure may need to be performed numerous times to find the optimum number of 
coefficients to meet the requirements of the application. 

For this simulation, Square Root Raised Cosine filters will be used on the 
transmitting and receiving end of the system. These filters are defined by the frequency 
response in Equation 3.3. Using FIR filter design, the filter can be designed directly from 
the filter response, whereas in UR design, an appropriate analog filter containing poles 
and zeros is needed to begin the design. In addition, Raised Cosine filters need linear 
phase to achieve zero ISI [12]. Using FIR filters, linear phase can be achieved. The only 
drawback to the design of Square Root Raised Cosine filter using FIR method is the large 
number of coefficients needed to achieve an accurate magnitude response. For this 
simulation FIR filter design method will be used to design the Square Root Raised Cosine 
filter. The filter coefficients were generated using a program written in Mathcad. 

3.3.3 FIR Filter Design of Transmit Filter 



For illustration, the design procedure for obtaining the FIR filter coefficients for 
the transmit filter will be demonstrated below. The procedure for obtaining the receive 
filter coefficients is exactly the same, the only difference being the filter transfer function. 

FIR filter design can either be done using the filter transfer function or the filter 
impulse response. The square root raised cosine filter is both specified in the frequency 
domain and in the time domain. The transmit filter needs to be cascaded with x/sin(x) in 
order to obtain zero ISI for a NRZ input, so it is more convenient to start in the frequency 
domain. The transfer function for the transmit filter is 


' K-f-T s 
sin(rc • f • T s ) 


for |f| < fj 


H(f) = 


rc f-T s , Jt-(f-fi) , 

cos ( — ) 

sin(7t • f • T s ) 4 • (B 0 - f t ) 


for fj < |f| < 2-B 0 - fj 


0 


for |f| > 2 • B 0 - fj 


(3.4) 


and is shown in Figure 3.3. The impulse response is obtained by taking the inverse fast 
Fourier transform of the filter transfer function. This analog transfer function first needs 
to be sampled. The signals in the simulation have a symbol rate R s , and symbol period 

of T s . In the time domain, the signal waveform is sampled at 16 samples/symbol (N ss = 
16). The sampling frequency is given by f s = R s • N ss Hz, and the spacing between 
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f 


Figure 3.3: Square root raised cosine filter with 
x/sin(x) equalization and rolloff factor a = 0.45 ( R s = l)(normalized frequency) 


samples is At = — — seconds. In the frequency domain, the filter response is periodic 

N ss 


in f s (see Figure 3.4). The total frequency span is f s • N ss , and for a FFT of length N, the 


f • N 

frequency spacing is Af = — — . After the filter transfer function has been sampled 

N 

with frequency spacing Af , the impulse response of the filter is obtained by taking the 
inverse FFT (IFFT). The impulse response of the square root raised cosine filter is shown 
in Figure 3.5. The impulse response of this filter is infinitely long. Truncation of the 
impulse response will allow us to have a finite number of coefficients in our FIR filter. 
The impulse response should be truncated at a point where the response has sufficiently 
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died down. In this case we will chose the number of one-sided coefficients (M) to be 


100 . 


H(f) 

o 

0 2 4 6 8 10 12 14 16 

f 

Figure 3.4: Transmit filter transfer function (normalized frequency) 




Figure 3.5: Filter impulse response 


The impulse response is actually two sided (see Figure 3.6). Because there is filter output 
before t = 0, the filter is non-causal. In order to correct this, a delay will be introduced to 
make it causal. This is done by shifting the impulse response M coefficients to the right 
(see Figure 3.7). In order to improve the design of this filter, windowing techniques will 
be used. Instead of abruptly cutting off the coefficients at ±M, window functions 
smoothly reduce the filter coefficients to zero. For this filter, a Hamming window is used 


65 




to accomplish this (See Figure 3.8). The filter coefficients are ‘windowed’ by the 
equation: 


h(n) = h ideal • w(n) n = 0,1 ... 2-M 
where the window function for a Hamming window is given by: 


w(n) = 0.54 - 0.46 • cos(— -) 

M 



Figure 3.6: Non-causal impulse response 



Figure 3.7: Causal impulse response 
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The resulting filter coefficients are used to filter a signal x(n) using the following 
equation: 

2*M 

y(n) = £h(k)-x(n-k) 
k = 0 

The output y(n) of the filter is dependent on the current input x(n), and the 2 • M previous 
inputs. There is no feedback involved, so the filter is always stable. The magnitude 
response of the filter as compared to the original filter is given in Figure 3.9. 



Figure 3.8: Hamming window function 


By using more filter coefficients, the response of the FIR filter will be closer to the 
original transfer function of the filter. 

3.4 Adding noise 

In order to make the simulation as realistic as possible, the channel model is 
chosen to be the Additive Gaussian White Noise (AWGN) channel. The noise signal can 
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be generated by producing a Gaussian, or normal sequence with standard deviation o 


using (3.2). 



f 


Figure 3.9: FIR frequency response vs. analog filter response 
for a square root Raised Cosine filter with x/sin(x) equalization 


3.4.1 Noise Equivalent Bandwidth 

The noise equivalent bandwidth B N will have to be calculated in order to add the 

correct amount of noise to the simulation. Consider the lowpass filter in Figure 3.10 a) 
with transfer function H(f). 
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Figure 3.10 (a) Receive filter (b) Noise equivalent bandwidth filter 


If white noise with a power spectral density 



is applied at the input to the filter, the 


total noise power at the output is 


P = 


00 M 

{^-•|H(f)| 2 df 


(3.4) 


= N 0 J|H(f)| 2 df ■ ■ (3.5) 

0 


3.4.2 Calculating the noise variance 

Now consider an ideal filter Hj(f) with single sided bandwidth B N as in Figure 

3.10 b). If the same white noise signal is applied to the input of the ideal filter, the total 
noise power at the output of the filter is 
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(3.6) 


XT 

Pideal = | -f !H[(f)| 2 df 

-B n 

= N 0 B n |H(0)| 2 (3.7) 


The bandwidth of the ideal filter can be chosen so that the total noise power at the output 
of the ideal filter is equal to the total noise power at the output of the real filter. Equating 
(3.5) and (3.7) and solving for B N yields 


J|H(f)| 2 df 

0 

|H(0)| 2 


(3.8) 


White noise has a constant PSD for all frequencies 

N n 

Snn( 0 = f° r - oo < f < oo 


S NN( f ) 


y 

X 

s. 

N 0 /2 

< -PK 

» 


Figure 3.1 1: Power spectral density of AWGN 


Unfortunately, this signal take an infinite amount of power to produce. For simulation 
purposes, we will be working with a system that has a finite bandwidth. The receive filter 
has a bandwidth B. The sampling frequency f s is chosen to be greater than 2B. If we use 

bandlimmited white Gaussian noise with a constant PSD over the simulation bandwidth, 
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< f < 


S N .N.(0 = ^ for 


the response of the system will be the same whether S NN (f) or S N N (f) is used [5]. 


H(f) 



a) 


^NsNsO ) 

A 

JL 

N 0 / 2 



b) 

Figure 3.12: a) Filter transfer function b) Bandlimited AWGN before filtering 


The total power for the bandlimited AWGN is equal to 

- 2 Npf s 

L) Y 


( 3 . 9 ) 


After AWGN is passed through a filter with noise equivalent bandwidth B N , the total 
noise power is equal to 


2 _ NO 


2-B n - N 0 B n - 


SNR 


( 3 . 10 ) 
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where P is the total signal power and SNR is the signal to noise ratio. Combining (3.9) 
and (3.10) and solving for g 2 


SNR 2 • Bi 


(3.11) 


G x 2 is the total noise power that is to be added to the channel before the Rx filter. The 
noise can be generated by using (3.2) to create a sequence of Gaussian random variables 
with jx = 0 and a = o x . 

Sometimes it is convenient to express the probability of a bit error P e as a 


function of p = where E b is the bit energy and N 0 is the noise density. The noise 
N 0 

variance will be expressed differently than (3.1 1). 

E b 

Let p = . In terms of N q 

N 0 


Nn = 


(3.12) 


The signal power P is equal to 


Eu -k 


(3.13) 


where k is the number of bits per symbol. From (3.13), the energy per bit is 


E b = 


(3.14) 


The total noise power before filtering is equal to g x 2 = — ^-5- from Figure 3 b). 


Combining this with (3.12) and (3. 14) yields 
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(3.15) 


a 


2 

x 


P T S f s 

2 • p • k 


The sampling frequency is defined as 


f s ~ P-s ‘ N ss 


so (3.16) can be rewritten as 


2 = P-N ss 

2 • p • k 


(3.16) 


(3.17) 


The noise can be generated by using (3.2) to create a sequence of Gaussian random 
variables with p = 0 and o = a x . 


3.4.3 Ideal Channel model 

An ideal channel can be used to obtain results that are not degraded by the filters. 
There is no need to simulate the analog signal, and modulate, filter, and add noise. The 
binary data can have Gaussian noise added directly to the data bits, and then these bits 
can either be soft or hard decoded. For each binary bit Uj , the noise corrupted bit can be 

obtained by using (3.2) with p = 1 for uj = 1, and p = -1 for Uj =0. The variance of the 

noise is 


o 


2 


1 

2p 


where p is the signal to noise ratio. For hard decisions Vj , if the noise corrupted Uj is 
greater than zero, then Vj = 1. Otherwise, Vj = 0. Using this method, ideal BPSK is 
obtained. 
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3.5 Simulation Results 


The simulation was run without coding systems to verify the simulation was 
working properly. The results were compared to the theoretical BPSK results, which are 
given by 


P p = — • erfc 
e 2 


W N o ) 


where E b is the energy per bit, N 0 is the single sided noise power spectral density, and 


erfc(x) = -T='J 


. -lT/2 


du. 


The results of the simulation are shown in Figure 3.13. Looking at the results of the 
simulation, two observations can be made. First, for greater values of E b /N 0 , the 

simulation results deviate from the theoretical curve. This can be attributed to the effects 
of filtering. Second, for higher values of rolloff, the performance of the simulation 
becomes closer to the theoretical. This can be attributed to the increase in bandwidth 
which comes from an increase in the rolloff value. 
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Chapter 4 


Erasure Methods and Simulation Results 

The purpose of this chapter is to investigate the performance of the concatenated 
system given in Figure 1.1, and to discuss and propose methods of declaring RS symbol 
erasures to improve performance. It is the goal of these methods to recognize symbols 
that are in error and erase them, thus utilizing the full capacity of the errors and erasures 
Reed Solomon decoder. The success of these methods is highly dependent on being able 
to successfully identify the errors in each Reed Solomon codeword (RSW). The 
performance of the concatenated system without erasures was tested. In addition, a basic 
method for erasure declaration using the reliability information provided by the SOVA 
was implemented. This method is based on [14]. Two new procedures for declaring 
erasures are proposed. Both use the reliability information from the SOVA in addition to 
information provided by successfully decoded RSW in the deinterleaving frame. The 
results of these two methods are presented and compared to Paaske’s [7] method. 

4.1 Concatenated system simulation results 

The performance of the concatenated system presented in Figure 1.1 was 
investigated. More specifically, the concatenated system used by NASA for deep space 
communication has been simulated. This system uses a rate 1/2, K = 7 convolutional 
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code with Viterbi decoding as the inner code and a (255, 223) Reed Solomon code as the 
outer code. This coding system was tested using a real system containing square root 
raised cosine FIR filters using a rolloff value of a = 0.5, and an ideal system. The 
performance evaluation studied the effects of using no interleaving, and interleaving 
depths I = 2, 6, 8. In addition, the effect of increasing the truncation lengths from 8 = 32 
to 8 = 100, and soft decision levels from L = 8 to L = 64 were simulated. 

Five systems have been developed to study the effects of varying the truncation 
length and number of soft decision levels, and are presented in Table 4.1. These systems 
will be simulated for no interleaving, and for interleaving depths 1 = 6 and 8. For the 

following results, the gain is measured at a bit error rate of 10 -5 . 


Table 4. 1 Various systems simulated 



Real/Ideal 

Truncation 
length 8 

Soft decision 
levels L 

System 1 

Real 

32 

8 

System 2 

Ideal 

32 

8 

System 3 

Ideal 

32 

64 

System 4 

Ideal 

100 

8 

System 5 

Ideal 

100 

64 


Figure 4. 1 contains the simulation results for the concatenated coding system with 
no interleaving used. It is shown that the real system with 8 = 32 and 8 level soft decision 
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Viterbi decoding (system 1) results in a loss of about 0.3 dB when compared with the 
ideal system with the same truncation length and soft decision levels (system 2). There is 
approximately a gain of 0.2 dB when the decision depth 5 is increased from 8 = 32 to 8 = 
100. This can be illustrated in Figure 4.1 by comparing the difference in performance 
between system 2 (8 = 32, L = 8) and system 4 (8 = 100, L = 8). Likewise, there is also a 
0.2 dB gain when the decision depth in system 3 is increased from 8 = 32 to 8 = 100 to 
form system 5. When the number of soft decision levels is increased from L = 8 in 
system 2 to L = 64 in system 4, the result is about a 0.15 dB gain. Likewise, the increase 
in soft decision levels from system 3 to system 5 results in roughly 0.15 dB gain. The 
combined use of 8 = 100 and L = 64 results in about 0.35 dB gain over the ideal system 
using 8 = 32 and L = 8. These gains are consistent for interleaving depths 1 = 6 (Figure 
4.2) and 1 = 8 (Figure 4.3). 

The benefits of using an interleaver is demonstrated in Figure 4.4, where the 
simulation results for system 5 using no interleaving, and depths I = 2, 6, and 8 are 
presented. It has been shown that the use of interleaving depth 1 = 2 can provide 
approximately 0.2 dB gain over the system using no interleaving. When the interleaving 
depth is increased from I = 2 to I = 6, there is an additional gain of slightly more than 0.2 
dB, for a total of roughly 0.4 dB over the system using no interleaving. The use of an 
interleaver with depth 1=8 provided a minor gain over the system using 1 = 6. The total 
gain of using 1 = 8 when compared to the non-interleaved system is approximately 0.45 
dB. It should also be noted that the concatenated coding system using interleaving depths 
greater than 8 resulted in little or no improvement over the system using depth 1 = 8. 
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As was expected, gains can be obtained in the concatenated coding system by 
increasing the truncation length and the number of soft decision levels used. The use of 
an interleaver between the convolutional code and the Reed Solomon code also provided 
gain over the system using no interleaving. The gains obtained through these 
improvements do not require any additional coding or bandwidth expansion, only 
additional hardware size and complexity. 

Another way gain can be obtained without the use of additional coding or 
bandwidth expansion is through the use of erasures. A Reed Solomon code with 
minimum distance d min can correct v errors and p erasures as long as the inequality 

d min < 2v + p is satisfied. Clearly, if errors are transformed into erasures, 
performance can be improved. Because the error positions are not known ahead of the 
decoding process, it is necessary to find methods that can identify unreliable Reed 
Solomon symbols. Once these symbols are identified, they can be erased. There is the 
possibility that correct symbols may be erased, however. For the performance to be 
improved, the erasure method must erase more errors than correct symbols. 

For the concatenated system in Figure 1.1, the Viterbi decoder produces hard 
output Uj 6 (0, 1}. It has been shown in section 2.3.3 that the soft output Viterbi 

algorithm produces a reliability value Lj associated with each outgoing bit Uj . This 

reliability information can give an indication to which Reed Solomon symbols are in 
error. In the first method investigated, bit reliabilities are transformed into RS symbol 
reliabilities. A table of the least reliable symbols for each RS codeword can then be 
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compiled, and the least reliable symbols can be erased. The first method for declaring 
erasures in presented in greater detail below. 

4.2 Erasure Method 1 

The first method for declaring RS symbol erasures is based on the method found 
in [14]. The SOVA is used to obtain the reliabilities of each output bit from the Viterbi 
decoder. The output of the SOVA is quantized using q+1 bits, where q bits represent the 
magnitude of the decision (the reliability), and 1 bit represents the hard decision. In 
general, these reliabilities are real valued, and have a range between 0 and infinity. It is 
not necessary to quantize the reliability values between 0 to infinity, but rather between 0 
and some value L max . The reliability values produced by the SOVA will be assigned one 

of 2 q levels between 0 and L max . If a reliability is greater than L max , then the 

reliability level is set to 2 q (i.e. the maximum level). For this simulation q = 8 and 
L ma x = 8 were used. From here on out, the “reliability value” Lj for a bit will refer to 

the q-bit quantized level. 

Before declaring RS symbol erasures, the bit reliabilities from the output of the 
SOVA need to be converted to RS symbol reliabilities. This can be accomplished by 
simply using the minimum bit reliability in a symbol as the symbol reliability. This can 
be rationalized as follows. If a bit has a small reliability, and is in error, the symbol will 
be in error also. Therefore, the best information contained in the symbol reliability is 
contained in the minimum bit reliability. For each RS codeword (RSW), a reliability 
table (RT) is formed. The RT has m e positions, where m e is the maximum number of 
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erasures allowed. If the least reliable symbols in the RSW occur at positions kj, where j 
= 1, 2, m e , then RT(j) = kj. RT(1) contains the least reliable position, RT(2) 

contains the 2nd least reliable symbol, etc. The RS decoder will first attempt to decode 
without declaring erasures. If this is a successful decoding, stop. If not, the two symbols 
with the smallest reliabilities are erased (RT(1) and RT(2)), and decoding is attempted 
again. If successful, decoding stops. If unsuccessful, the four least likely symbols are 
erased (RT(1) through RT(4)) and decoding is attempted again. For each unsuccessful 
decoding, two more erasures are declared until either a successful decoding is achieved, 
or the maximum number of erasures m e is reached. This maximum number is chosen to 

be 16 erasures. Using more than 16 erasures gave poorer results due to the fact that the 
symbols reliabilities at the output of the SOVA hit more correct symbols than errors. In 
addition, the probability of a decoding error is also increased. A decoding error occurs 
when a codeword contains more than t errors, and fails to notice that it does. The decoder 
claims that there are less than t errors. This is different than a decoding failure. A 
decoding failure happens when there are more than t errors, and the decoder detects that it 
contains more than t errors. For a t error correcting RS code, the probability of a 
decoding error is 1/t! [7, 14, 19]. The procedure for erasure decoding using Method 1 is 
presented below. 

4.2.1 Procedure for erasure decoding using Method 1 

The decoding of each RS codeword involves the following steps: 

1. Set the number of erasures n e =0. 
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2. Attempt to decode the RS codeword using errors only decoding. If successful, stop. 

3. Set n e = n e + 2. 

4. Erase the n e least reliable symbols as determined by the reliability table, and attempt 
error and erasures decoding. If successful, stop. 

5. If n e = m e (the maximum number of erasures allowed), then stop. 

6. Go to step 3. 

7. Stop. 

4.2.2 Simulation Results for Method 1 

The results of the simulation for Method 1 are presented in Figures 4.5 through 
4.7. These figures contain the results using system 1, 2, and 5 for no interleaving (Figure 
4.5), and interleaving depths I = 6 (Figure 4.6) and 1 = 8 (Figure 4.7). The results of the 
corresponding concatenated code using no erasures is also presented for comparison 
purposes. From these figures, it is evident that the use of erasure Method 1 results in 
approximately 0. 1 dB gain, over the non-erasure concatenated system. This result was 
consistent regardless of the truncation length, number of soft decision levels, and 
interleaving depth used by the concatenated system. It should be noted that the gains 
obtained by using different interleaving depths with erasure Method 1 were identical to 
the gains obtained when using different interleaving depths in the standard concatenated 
system. In addition, it should also be noted that increasing the decision depth and number 
of soft decision levels for the systems simulated using erasure Method 1 were consistent 
with the results for the concatenated system presented in section 4.1. 
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Figure 4.7 Simulation results for Method 1 with interleaver depth 1 = 8 


It has been shown that the use of erasure Method 1 results in approximately 0. 1 
dB gain over the concatenated system using no erasures. The reason that the performance 
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is not greatly increased is because the erasure declaring process is not optimal. In general, 
when erasing two symbols at a time, gain is only achieved when both erased symbols are 
errors (GE). If only 1 erasure is a GE, the capacity of the code stays the same, and if 
neither erasure is a GE, then capacity of the code actually decreases. Method 1 uses the 
reliability information generated from the SOVA, and the symbols with the smallest 
reliabilities are systematically erased. A symbol with a small reliability does not 
guarantee a GE. For example, the smallest reliability value ( Lj =0) occurs when there is 

a metric difference A = 0. Because the path metrics are equally likely, the probability the 
wrong path has been chosen is 0.5. This is why a small reliability value does not 
guarantee a wrong path has been taken. The reliability table used in Method 1 can be 
modified to give more reliable information. This can be accomplished by using decoded 
Reed Solomon codewords in the deinterleaving frame to provide information to the non- 
decoded Reed Solomon codewords. This information can be used to construct a 
reliability table with more accurate information, and then the unreliable symbols can be 
converted into RS symbol erasures. The second method for declaring erasures is 
presented below. 

4.3 Erasure Method 2 

It has been found that the burst errors at the output of the SOVA contain the same 
reliabilities within the burst. Consider the two paths merging in state S k at time t in 

Figure 4.8. The survivor path has metric Mj and the competing path has metric M 2 . 
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Figure 4.8 The updating of symbol reliabilities in the SOYA 


The reliabilities Lj on the surviving path in the SOVA are updated by selecting the 
minimum between Lj and A, where A is the difference between the two paths merged in 
state S k . When the bit reliabilities are converted to symbol reliabilities, the minimum bit 
reliability is chosen as the symbol reliability. Assuming that A is the minimum and this 
is part of the surviving path, then consecutive symbols in the surviving path have 
identical reliabilities. In the standard Viterbi decoder, error events occur when at time t 
the decoder chooses the wrong path. This is also what happens in the SOVA. If the path 
with the maximum metric is the wrong path, then all of the reliabilities for the differing 
bits will be equal to A. If this is the path selected by the SOVA as the output of the 
decoder, then the reliability will be equal for the length of the burst. It is apparent from 
Figure 4.8 that all symbol reliabilities in the erroneous path are equivalent. This 
information can be used to distinguish error paths and declare erasures, and can be 


91 



accomplished by using a method similar to the one proposed by Paaske. Assume that 
after an initial decoding attempt of a deinterleaving frame, some of the RSW have been 
successfully decoded, while others have not. The error positions in the successfully 
decoded codewords are known. Using the above observation, if a neighboring symbol in 
an undecoded RSW has an identical reliability value as the corrected symbol in a decoded 
RSW, it is highly probable that that symbol is also an error and can be erased. 

As mentioned previously, a symbol with a low reliability value does not 
necessarily denote a symbol error. From Figure 4.8, it can be seen that if A is small, then 
it may be contained in the reliability table for a non-decoded RSW. Assume that a RSW 
is decoded and the error positions are known. The correct positions are also known. 
These correct positions can be used to eliminate other correct symbols from being 
included in the reliability table for non-decoded RSW. If the neighboring positions of a 
correct symbol have the same reliability, then it is highly probable that this symbol is 
correct also. The procedure for erasure decoding using Method 2 is summarized below. 

A deinterleaving frame with interleaving depth I contains I RSW. Decoding of 
each RSW in the frame is first attempted using Method 1 . For each RSW, a table of the 
least reliable symbols is compiled, and for each unsuccessful decoding attempt, two 
additional symbols are erased, and decoding is attempted again. This process continues 
until either a successful decoding occurs, or a maximum amount of erasures has been 
reached. If less than I codewords are successfully decoded, redecoding is attempted. The 
RT used in Method 1 contains the least reliable symbols as determined by the output of 
the SOYA. The reliability table can be modified using information provided by the 
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decoded RSW in the deinterleaving frame. Let RSW k (i) denote the k-th symbol of the i- 
th RSW in the deinterleaving frame, and let L k (i) be the reliability associated with 
RSW k (i) as determined by the SOVA. For convenience, if i + j > I, then RSW k (i + j) 
corresponds to RSW k+1 (i + j - 1) and if i + j < 1, RSW k (i-j) corresponds to 
RSW k-1 (i - j + 1) . To simplify notation, RSW k (i ± j) will be used even if the position 
is at k + 1 . In addition, each symbol will have a flag associated with it. Let F k (i) be the 
flag for the k-th symbol in the i-th RSW. There are three possible values for F k (i) : 

F k (i) =0 Unknown 

F k (i) = 1 Possible Good Erasure (PGE) 

F k (i) = 2 Possible Bad Erasure (PBE) 

Using the already decoded RSW in the deinterleaving frame, the flags for all the symbols 
in the undecoded RSW will be updated. The table of least reliable symbols can be 
modified using the symbols that are flagged as a PGE (F k (i) = 1) as the least reliable. 

The remainder of the table is filled with the symbols flagged as unknown (F k (i) = 0) 

with the minimum reliabilities. Symbols that are presumed to be correct are flagged as a 
PBE ( F k (i) = 2) to avoid being used in the RT. It was found through simulation that the 

probability that a symbol flagged as a PGE is a GE is 0.93 and the probability that a 
symbol flagged as a PBE is a BE is 0.99. 

Figure 4.9 shows a partial deinterleaving frame, two bursts, and the reliability 
values associated with each symbol. Note that the reliability in each of the burst errors is 
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identical. Assume that RSW(i) has been successfully decoded in the deinterleaving frame 
of Figure 4.9. After the decoding, the error positions in RSW(i) are known. The error at 
position k has a reliability of L k (i) = 2. Erasure Method 2 checks the reliability at 
position k in RSW(i-l). If the two have equal reliabilities (i.e. L k (i - 1) = L k (i)), then 
RSW(i-l) position k is flagged as a PGE ( F k (i — 1) = 1). If the two are not equal, 

searching in this direction stops, and searching in the other direction begins. Position k in 
RSW(i+l) is checked for equal reliabilities, and if so, position k in RSW(i+l) is flagged a 
PGE. Searching RSW(i±j) continues in both directions until either L k (i±j) * L k (i), 

or j =1-1. The later condition can be reasoned by noting that position k in RSW(i+I) is 
simply position k+1 in RSW(i). 



Figure 4.9 Deinterleaving frame with reliability information 


Assume that the two least reliable symbols contained in the RT for the first 
decoding try (Method 1) are symbol k with reliability L k = 2 and symbol k + 2 with 
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reliability L k+2 = 11. If two erasures are declared, position k would be a GE and 
position k + 2 would be a BE. After decoding RSW(i), it has been found that position k + 
2 is correct. Because L k+2 is a small reliability, it is most likely contained in the RT for 

each RSW that has L k+2 = 11. The same method used to flag PGE can be used to flag 

PBE. After doing so, position k+2 will not be in the RT for the next iteration of 
decoding. The procedure used to update the flags (UFP) in the undecoded RSW using 
information in decoded word RSW(i) is as follows: 

4.3.1 Procedure for updating the flags (UFP) 

1. Set k = 0 

2. If RSW k (i) is an error set F k (i) = 1 (PGE). Else, set F k (i) = 2 (PBE). 

3. If L k (i) > L max go to step 13 

3. Initialize j =1. 

Start looping backwards. RSW k (i - j) j = 1 , 2, ..., I - 1 . 

4. If L k (i - j) L k (i) , go to step 8. 

5. If L k (i - j) = L k (i) set 

F k (i-j) = F k (i) 

j=j + l 

6. If j = I - 1, go to step 8. 

7. Go to step 4. 

Start looping forwards. RSW k (i + j) j = 1, 2, ..., I - 1. 
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8. Initialize j = 1. 

9. If L k (i + j) * L k (i) , go to step 13. 

10. If L k (i + j) = L k (i) set 

Fk(‘ + j) = F k(0 
j=j + l 

11. If j = I - 1 , go to step 13. 

12. Go to step 9. 

13. k = k+ 1. 

14 If k = 254, Stop else go to step 2. 

After the flags have been updated using the information provided by the 
successfully decoded RSW, the RT for each of the yet to be decoded RSW needs to be 
updated. The procedure for updating the RT (URTP) is presented below. 

4.3.2 Procedure for updating the reliability table (URTP): 

1. Determine the number of PGE ( n PGE ) in RSW(i) by checking the symbol flag F k (i) . 
These GE occur in positions k: , j = 1,2,..., n PGE . 

2. Fill the first n GE positions of the reliability table with the positions where a GE has 
been flagged. RT(j) = kj for j = 1,2, n GE 

3. If n PGE < m e , fill the remainder of the RT with the symbols with the minimum 
reliabilities determined by the SOVA output. Out of all the symbols that have F k (i) = 0, 
determine the m e - n GE minimum symbol reliabilities. These minimum values occur in 
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positions k^, t - 1, 2, m e - n PGE . Fill the remaining m e - n PGE positions in RT 
with the symbols with the minimum reliabilities. RT(n GE + 1) = k^, i - 1, 2, m e - 

n PGE * 

The procedure for decoding using Method 2 is presented below. 

4.3.3 Procedure for decoding using Method 2 

1. Attempt initial decoding using Method 1. 

2. Set i d = number of correctly decoded RSW. 

3. If i d = I, go to step 9. 

4. If i d = 0, go to step 9. 

5. For each decoded codeword, declare flags in the undecoded RSW using UFP. 

6. For each non decoded codeword, update the reliability table using URTP. 

7. Attempt decoding using Method 1 with the updated RT. 

8. If step 7 yields at least one successfully decoded RSW, go to step 2. 

9. Stop. 

4.3.4 An example of erasure decoding using Erasure Method 2 

As an example, consider the deinterleaving frame in Figure 4.10. This frame was 
simulated at a signal to noise value of E b / N 0 = 1.9 dB. The error positions are denoted 

by the shaded areas, and the reliability values are given. This frame contains 8 RSW with 
a total of 161 symbol errors. It should be noted that each RSW in the deinterleaving 
frame contains more than 16 errors. If no erasure information was used, every RSW in 
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the frame would fail to decode. The reliability table for each RSW in the frame for the 
first iteration is presented in Figure 4.1 1, where the reliability value and symbol position 
are given (L k , k) . In addition, the number of decoding trials for each codeword is also 

shown. The RT is initially formed by using the minimum reliabilities generated by the 
SOVA. For the first iteration, decoding of each RSW is attempted using Method 1. 
RSW(l) successfully decodes after 6 decoding trials using 10 erasures. RSW(3) is 
successfully decoded after declaring 4 erasures (3 decoding trials), and RSW(8) decodes 
using 2 erasures (2 decoding trials). All other RSW fail to decode after 9 decoding trials 
apiece. Now, the RT is modified by updating the flags using UFP, and the updated RT is 
compiled using URTP giving the symbols with flag = 1 (PBE) the highest priority in the 
table. The modified RT is in Figure 4.12. Decoding is attempted again using Method 1 
with the updated RT. For this iteration, RSW(2) decoded successfully using 14 erasures 
(8 decoding trials), RSW(4) using 12 erasures (7 trials), and RSW(7) using 8 erasures (5 
trials). Decoding attempts for RSW(5) and RSW(6) are once again unsuccessful. The 
flags for RSW(5) and RSW(6) are updated using the information provided by RSW(2), 
RSW(4), and RSW(7). The updated RT is generated using the updated flags, and is 
presented in Figure 4.13. Decoding is attempted again using Method 1. RSW(5) is 
decoded using 16 erasures (9 decoding trials), and RSW(6) is decoded using 14 erasures 
(8 decoding trials). Erasure Method 2 required a total of 111 decoding trials to 
successfully decode this frame. This results in an average of 13.9 decoding trials per 
RSW. 
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Figure 4.10 An example of a deinterleaving frame with reliability values 
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Figure 4. 1 1 The reliability table for the first iteration 
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Figure 4. 12 The reliability table for the second iteration 
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Indicates an erasure hits an error (GE) 


Figure 4. 13 The reliability table for the third iteration 


4.3.5 Results for Erasure Method 2 

The statistics for Method 2 are presented in Table 4.2. The RS symbol error rate 
at the output of the Viterbi decoder is given. The percentage of frame failures and RSW 
failures after errors only RS decoding gives an indication how the standard concatenated 
system performs. The percentage of frames failures and the percentage of RSW failures 
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using Method 2 are also compiled for various values of E b /N 0 . The number of 
decoding trials per RSW is also given. 


Table 4.2 Simulation results using Method 2 with interleaving depth 1 = 8 



Without Erasures 

Using Method 2 


Byte Rate 
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% Frames 
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E b /N 0 ’ 
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When the results in Table 4.2 are compared to Paaske’s results, the first difference 
that is noticed is the difference in the RS symbol error rate at the output of the Viterbi 
decoder. The symbol error rate given in Table 4.2 is obtained using system 5 (5 = 100 

and L = 64). This symbol error rate is approximately 3 • 1CT 3 lower than the symbol error 
rate in Paaske’s results. The better performance at the output of the Viterbi decoder is 
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most likely due to the increased soft decision levels (L = 64) as compared to the 8 level 
soft decision used by Paaske. This slightly smaller RS symbol error rate translates into 
slightly better results for the errors only decoding. The frame error rate and RSW error 
rate in Table 4.2 are typically 2% or 3% smaller than the results that Paaske obtained. 

The results obtained using Method 2 are better than Paaske’ s results in three 
regards. The percentage of RSWs and frames in error is lower. For example, at E b / N 0 

= 1.9 dB, the percentage of frames and RSWs in error using Paaske’s method are 26.8% 
and 16.3%, compared to 20.2% and 8.3% using Method 2. Method 2 also obtains these 
results using considerably less decoding trials. Paaske’s method requires an average of 
1 12.4 decoding trials per RSW to obtain the reduction. Method 2 requires an average of 
7.7 decoding trials per RSW. The third improvement is in the ability to obtain 
improvement for lower values of E b /N 0 . The results in [7] indicate that no 

improvement is obtained at E b / N 0 = 1.8 dB. Using Method 2, the percentage of frames 
in error was reduced from 99.4% to 61.1% and the percentage of RSW in error was 
reduced from 79.9% to 36.6%. The average number of decoding trials per RSW needed 
to accomplish this is 13.2. Error reductions were obtained for values of E b / N 0 as low 

as 1.7 dB. Both Method 2 and Paaske’s method obtained 0.3 dB gain over the system 
using no erasures. 

The reduction in average number of RSW trials as compared to Paaske’s method 
can be attributed to the ability to use erasures in the first decoding attempt. This is useful 
for decoding RSW that contain more than 16 errors, as in the example in Figure 4.10. 
Each RSW in the frame contains more than 16 errors. Paaske’s method would obtain no 
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successful decoding in the first iteration, and would have to resort to using EP4 to obtain 
the first successful decoding. EP4 randomly erases two symbols and attempts decoding. 
This repeats until a successful decoding or this has been attempted T max times. It can 

easily be seen that this frame would require a large amount of decoding trials. In 
addition, Paaske’s method requires two successfully decoded RSW with errors in the 
identical positions in order to use EP1. Erasures declared using EP1 have a probability of 
0.96 of being a GE. If not, a less reliable procedure must be used (EP2-EP4). Erasure 
Method 2 only requires 1 RSW to be able to declare erasures with a probability of 0.93 of 
being a GE. It can be seen that erasure Method 2 converges on the erasures quicker than 
Paaske’s method, and thus, requires less decoding trials per RSW. 

4.4 Erasure Method 3 

Method 2 can be modified to reduce the number of decoding trials per RSW by 
making a few observations. In Method 1, erasures were erased two at a time until either a 
decoding success was obtained, or a maximum number of erasures were declared. This 
was done because symbols with low reliabilities were not guaranteed to be GE. Erasing 
more than two at a time may cause the decoding capability of the code to be decreased. 
This occurs if the number of BE is greater than the number of GE. In Method 2, when a 
PGE is declared, it is highly probable that this is a GE (0.93). Instead of erasing two 
symbols from the RT at a time, all symbols flagged as PGE are automatically erased. If, 
after decoding with these symbols erased it is still not successful, then two more symbols 
are erased until the maximum allowable erasures is reached. 
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A second modification updates the flags after a RSW has been decoded correctly, 
rather than after all of the RSW in the frame have been attempted. This allows for the 
possibility for the highly reliable PGE to be declared in the first decoding pass. Both of 
these modifications reduce the average number of RS decoding trials when compared to 
Method 2, with very little effect on the performance. The procedure can best be 
demonstrated with an example. 

4.4.1 An example of erasure decoding using Erasure Method 3 

Consider the deinterleaving frame in Figure 4.10. The reliability table before the 
first decoding attempt is given in Figure 4.1 1. Much like Method 2, decoding of RSW(l) 
is first attempted using erasure Method 1 (i.e. decoding is attempted using no erasures, 
and decoding is repeated until either a successful decoding or a maximum number of 
erasures has been reached). RSW(l) is finally decoded after 6 decoding trials and 10 
erasures. In Method 2, the next step would be to attempt decoding of RSW(2) using 
Method 1. Method 3 instead updates the flags immediately after RSW(l) has been 
successfully decoded. This results in 13 PGE being declared in RSW(2), 12 of which are 
GE. This modified reliability table is given in Figure 4.14. At this point, RSW(2) is 
ready to be decoded. Instead of erasing two symbols at a time as with Method 2, all PGE 
are erased and decoding is attempted. For RSW(2), this still does not yield a successful 
decoding, but as can be seen from Figure 4.14, erasing two at a time would achieve the 
same results, but with more decoding trials. Because the number of PGE are less than the 
maximum number of erasures allowed, two more symbols are erased, and decoding is 
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attempted. RSW(2) finally decodes once 16 erasures have been declared. This is 
accomplished in 3 decoding trials, as compared the 17 trials needed to decode RSW(2) 
using Method 2. The flags are now updated using the information provided by RSW(2), 
and the resulting reliability table is given in Figure 4.15. The flags declared from 
RSW(l) and RSW(2) yield 10 PGEs in RSW(3), 8 of which are GE. It should also be 
noted that the PBE that were declared removed potential BE from the table. Symbols at 
positions 207 and 112 would have been included in the reliability table. These would 
have resulted in BE if used. The decoding of RSW(3) is accomplished by erasing the 10 
PGE in the table. This resulted in a decoding success in only 1 decoding trial. The flags 
are modified using the information provided by RSW(3), and the resulting reliability table 
is given in Figure 4.16. There are 11 PGE in RSW(4), and all are GE. This RSW 
requires all 16 erasures for a successful decoding, and is accomplished with 4 decoding 
trials. The flags are updated and the resulting reliability table is shown in Figure 4.17. 
There are 15 PGE in RSW(5) and 14 are GE. The codeword contains 23 errors, and 
cannot be decoded with the 16 erasures in the table. This codeword will require 15 out of 
16 erasures to be GE for a decoding success. Two decoding trials were attempted on this 
codeword. Decoding of RSW(6) is attempted by declaring erasures at the positions where 
the 9 PGE have been flagged. This does not yield a successful decoding, so erasures are 
declared 2 at a time until the maximum is reached, at which point a decoding success has 
not been obtained. This required 5 decoding trials. RSW(7) is successfully decoded by 
erasing all 12 positions with PGE flags, and decoding. RSW(7) required only 1 decoding 
trial to decode successfully. The flags are updated, and the resulting reliability table is 
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given in Figure 4.18. RSW(8) contains 16 PGE, and 13 are GE. There are only 17 errors 
in RSW(7), and is easily decoded in one trial by declaring 16 erasures. The flags are once 
again updated and the resulting reliability table is given in Figure 4.19. Decoding of 
RSW(5) is attempted a second time. New information has been provided by RSW(7) and 
RSW(8). There are 17 PGE in RSW(5), and all are erased in the first decoding attempt. 
In Method 3, more than 16 erasures can be declared if all are PGE. RSW(5) successfully 
decodes and the updated reliability table is given in Figure 4.20. RSW(6) has 20 PGE, 
with 19 of these being GE. All 20 symbols are erased, and the word is successfully 
decoded. The total number of trials for this frame is 25, which corresponds to an average 
number of trials per RSW of 3.13. This is significantly less than the number of trials 
required for erasure Method 2 to decode the same frame. The results of the simulation are 
given in Table 4.3. It can be seen from the table that a gain of 0.25 dB is obtainable. 
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Figure 4.14 Reliability table after RSW(l) successfully decodes 








































































































Indicates an erasure hits an error (GE) 


Figure 4. 15 Reliability table after RSW(2) successfully decodes 
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Indicates an erasure hits an error (GE) 


Figure 4.16 Reliability table after RSW(3) successfully decodes 
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Indicates an erasure hits an error (GE) 


Figure 4. 19 Reliability table after RSW(8) successfully decodes 
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Figure 4.20 Reliability table after RSW(5) successfully decodes 
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Table 4.3 Simulation results using Method 3 with interleaving depth 1 = 8 



Without Erasures 

Using Method 3 


Byte Rate 

% Frames 

% RSW 

% Frames 

% RSW 

# of trials 

E b /N 0 

after VD 

in error 

in error 

in error 

in error 

per RSW 

1.7 

0.093 

100.0 

95.2 

94.7 

78.1 

9.7 

1.75 

0.086 

99.8 

89.7 

83.1 

59.2 

8.56 

1.8 

0.079 

99.4 

79.9 

64.1 

37.7 

6.68 

1.85 

0.072 

97.0 

66.1 

42.8 

20.5 

4.54 

1.9 

0.066 

91.6 

52.3 

23.6 

9.6 

3.03 

1.95 

0.061 

81.6 

37.3 

13.2 

4.6 

2.19 

2.0 

0.055 

66.8 

24.9 

6.7 

1.9 

1.51 

2.05 

0.05 

46.9 

14.9 

2.8 

0.6 

1.27 

2.1 

0.046 

32.7 

8.5 

0.6 

0.08 

1.13 

2.15 

0.042 

19.4 

4.4 

0.2 

0.0 

1.07 

2.2 

0.038 

10.1 

1.9 

0.0 

0.0 

1.03 

2.3 

0.031 

1.0 

0.2 

0.0 

0.0 

1.0 

2.4 

0.025 

0.2 

0.0 

0.0 

0.0 

1.0 


As can bee seen comparing the results of the two methods in Tables 4.2 and 4.3, 
the significant difference is the reduction in the average number of decoding trials needed 
per RSW when using Method 3. The tradeoff in this reduction in the average number of 
decoding trials is a slight increase in the percentage of frame and word errors using 
Method 3. Method 2 outperforms Method 3 by a few percent. The BER curves for 
Method 2 and Method 3 are given in Figure 4.21. Using erasure Method 2 and Method 3 
resulted in approximately 0.3 dB and 0.25 dB gain respectively over the concatenated 
system using no erasures. 
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Figure 4.21 Simulation results of the concatenated system using 
various erasure methods 
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Chapter 5 


Conclusions 

Performance of the concatenated coding system through the use of RS symbol 
erasures has been demonstrated. The first method investigated uses reliability 
information generated by a modified Viterbi decoder. This information is derived from 
the metric difference of two paths merging in each state. The reliability for the bits along 
the surviving path are updated using this metric difference. This method yielded about a 
0.1 dB improvement over the concatenated system using no erasures. This gain was 
independent of the truncation length, interleaving depth, and number of soft decision 
levels. 

In the second method proposed, the reliability table is refined using information 
provided by the decoded RSW in the deinterleaving frame. The error positions are 
known in the decoded word, and the method searches for equal reliabilities in 
neighboring, non-decoded RSW. If the reliabilities are the same, then symbols are 
erased. This method yielded approximately 0.3 dB gain over the standard concatenated 
system. The average number of decoding trials per RSW is substantially less than the 
results presented by Paaske [7]. Using Method 3, the number of trials per RSW is 
reduced even further. In addition, Methods 2 and 3 perform reasonably well at very low 
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values of Eb/No (i.e. < 1.9 dB). The cost is the increased complexity of having to use the 
SOVA in place of the standard Viterbi decoder. 

It has been shown that the use of convolutional code as the inner code and a Reed 
Solomon code as the outer code provide considerable gains. The rate 1/2 K = 7 
convolutional code concatenated with a (255, 223) t = 16 error correcting Reed Solomon 
code yielded about 8 dB gain over uncoded BPSK and approximately 9 dB with an 
interleaving depth of I = 8. The use of an interleaver with depth 1 = 8 gives a gain of 
approximately 0.5 dB as compared to the same system using no interleaver. It has been 
demonstrated that the use of a real system results in approximately a 0.3 dB loss when 
compared to the ideal system using the same truncation length and number of soft 
decision levels. It has been shown that a 0.15 dB improvement can be obtained by 
increasing the number of soft decisions used from L = 8 to L = 64. A gain of 0.2 dB is 
attainable if the truncation length in the Viterbi decoder is increased from 5 = 32 to 5 = 
100. These gains were obtained regardless of the interleaving depth used. 

5.1 Future Research 

Erasure method 1 declares RS symbol erasures by using reliability information 
provided by the SOVA. This method could be improved by using the maximum a 
posterori (MAP) algorithm in place of the SOVA. The MAP algorithm is an alternate 
method for decoding convolutional codes. The Viterbi decoder finds the maximum a 
posterori probability for the entire path through the trellis. The MAP algorithm, on the 
other hand, finds the maximum a posterori probability for each outgoing bit. The MAP 
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Algorithm was designed to minimize the word error probability rather than the sequence 
error probability as the Viterbi algorithm does. This algorithm is more complex than the 
Viterbi, but easily provides log-likelihood (reliability) values at its output. The reliability 
values generated by the MAP algorithm yield superior results when compared to the 
SOVA, but is very complex. The Map algorithm could be used in place of the SOVA in 
the concatenated coding system presented here. 

The performance of Method 2 could also be improved upon. The special case 
where GE are flagged BE, and vise versa, might be able to be detected and solved. 
Several methods were tried unsuccessfully to fix the problem. Performance 
improvements could be obtained by incorporating some of the elements used in Paaske’s 
method, such as EP1, EP3, and EP4. In Paaske’s method, when the highly reliable 
methods (EP1 and EP2) can not decode some of the RSW in the frame, the backup 
methods are used to exhaustively try and decode, which results in a high average number 
of decoding trials per RSW. But if through the use of this exhaustive search, a previously 
undecodable RSW is now correctly decoded, then this RSW can possibly help in the 
decoding of other non-decoded RSW in the frame. This is where a large number of 
decoding trials may be justifiable. When Methods 2 and 3 fail to decode a frame, there is 
no backup method that might provide improvement. The lack of a backup procedure, is 
one of the reasons that the number of trials is low compared to Paaske’s method. This 
lack of a backup procedure is somewhat justified. In Paaske’s method, EP4 is used to 
systematically guess the positions of the errors. This guessing is effective when the 
number of errors is not much larger than 18 or 19. For Methods 1 and 2, most RSW with 
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18 or 19 errors can be decoded using the initial decoding using Method 1. The SOVA 
output provides the decoder with enough useful information to decode some of the RSW 
that Paaske’s method would require EP4 to decode. For example, in section 4.5, there is 
an example RS frame decoded using Method 2. Every RSW in the frame contains more 
than 16 errors. If Paaske’s method was used to decode this frame, then EP4 would have 
to be used to obtain the first decoded RSW. Using Method 2, three RSW were 
successfully decoded on the first trial, and the others were successfully decoded using 
information provided by these decoded RSW. 

Another possible improvement is the use of iterative decoding. Iterative decoding 
takes the output of the Reed Solomon decoder and feeds this to the input to the Viterbi 
decoder. Paaske [7] used iterative decoding in his paper. If there are still undecodable 
RSW in the deinterleaving frame, the corrected codewords are sent back to the input to 
the Viterbi decoder. The corrected bits give the Viterbi decoder some of the states that 
the correct path took. Forcing the Viterbi decoder through these known states may help 
the Viterbi decoder in a better estimate of the correct path. The error rate at the output 
will be reduced, which will result in more RSW being decoded on the second trial. 
Iterative decoding could be used with Methods 2 and 3 presented in this report. There is 
the potential for larger gains using this repeated Viterbi decoding because the first 
iteration will have corrected some of the RSW, and provide the second Viterbi trial with 
additional information. It should be noted that this repeated Viterbi decoding is effective 
only when the output of the RS decoder contains some correctly decoded RSW. Methods 
2 and 3 allow for a reduction in RSW failure for Eb/No values at as low as 1.7 dB SNR. 
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Even though the reduction at the very low values of Eb/No is not very large, this little bit 
of reduction would help greatly if used in repeated Viterbi decoding trials. 

The systems investigated in this report could be simulated using SPW. SPW is a 
useful tool for modeling of communication systems. Many signal processing blocks such 
as filters, modulators, and channel models are contained in SPW’s library for use in the 
design of larger, complex systems. The SOVA and the errors and erasures RS decoder 
are currently not in any of the libraries found in SPW. SPW does allow for “custom 
coded blocks”. This tool can take C code, and create an SPW block based on this code. 
This custom coded SPW block can then be used in combination with other SPW blocks 
in the design of communication systems. 

The SOVA could be implemented in a single chip design. The method presented 
by Hagenauer and Hoeher could be used, or another method presented in [2]. 
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Appendix A 


Simulation Flow Charts 



Figure A. 1 Flow chart for erasure Method 1 
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Figure A.2 Flow chart for erasure Method 2 
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Figure A.3 Flow chart for the real simulation 
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Figure A.4 Flow chart for the ideal simulation 
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Appendix B 


Program Listing 


B.l Memory allocation functions 

int *ivector(long nh) 

{ 

int *v; 

v = (int *)calloc(nh, sizeof(int)); 
return v; 

} 


void free_2d_int_matrix(int row, int **a) 

{ 

int i; 

for(i=0;i<row;i++) 

free(*a); 

free(a); 

} 


int **int_matrix_2d(int row, int col) 

{ 

int i; 
int **a; 

a = (int **)calloc(row, sizeof(int *)); 
for(i=0;i<row;i++) { 

a[i] = (int *)calloc(col, sizeof(int)); 

} 

return a; 

} 


double **double_matrix_2d(int row, int col) 

{ 

int i; 

double **a; 

a = (double **)cal!oc(row, sizeof(double *)); 
for(i=0;i<row;i++) { 

a[i] = (double *)calloc(col, sizeof(double)); 

} 

return a; 

} 

void free_2d_double_matrix(int row, double **a) 

{ 
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int i; 

for(i=0;i<row;i++) 

free(*a); 

free(a); 

} 


double *dvector(long nh) 

{ 

double *v; 

v = (double *)calloc(nh, sizeof(double)); 
return v; 

} 


B.2 Random number generators 

double ran2(long *idum) 

/* Ran2 is a long period randum number generator ( > 2* 10 A 18). Returns */ 
/* a uniform random deviate between 0.0 and 1.0. Call with idum a */ 

/* negative integer to initialize; thereafter, do idum between successive */ 

/* deviates in a sequence. This subroutine is taken from the book */ 

/* "Numerical recipies in C" by Saul A. Teukolsky, William T. Vetterling, */ 
/* and Brian P. Flannery. */ 

{ 

int j; 
long k; 

static long idum2 = 123456789; 
static long iy = 0; 
static long iv[NTAB]; 
double temp; 

if (*idum <= 0) 

{ 

if (-(*idum) < 1) *idum=l; 
else *idum = -(*idum); 
idum2 = (*idum); 
for (j=NTAB+7; j>=0; j— ) 

{ 

k = (*idum)/IQl; 

*idum=IAl*(*idum - k*IQl) - k*IRl; 
if (*idum < 0) *idum += IM 1 ; 
if (j < NTAB) iv[j] = *idum; 


iy = iv[0] ; 

} 

k - (*idum)/IQl; 

*idum = IAl*(*idum - k*IQl) - k*IRl; 
if (*idum < 0) *idum += IM1; 
k = idum2/IQ2; 
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idum2 = IA2*(idum2 - k*IQ2) - k*IR2; 

if (idum2 < 0) idum2 += IM2; 

j=iy/NDIV; 

iy = iv[j] - idum2; 

ivfj] = *idum; 

if (iy < 1) iy += IMM1; 

if ((temp = AM*iy) > RNMX) return RNMX; 

else return temp; 

} 


float gasdev2(long *idum) 

/* This function returns a normally distributed deviate with zero mean and unit variance */ 
/* This function is more computationally efficient than gasdev 1 because there are no trig */ 
/* function calls and the funtion saves the extra deviate for the next funtion call. This */ 

/* subroutine is taken from the book ''Numerical recipies in C" by Saul A. Teukolsky, */ 

/* William T. Vetterling, and Brian P. Flannery. */ 

{ 

static int iset = 0; 
static float gset; 
float fac, rsq, vl , v2; 

if(iset == 0) 

{ 

do 

{ 

vl = 2.0*ran2(idum) - 1.0; 
v2 = 2.0*ran2(idum) - 1.0; 
rsq = vl*vl+v2*v2; 

} while (rsq >= 1 .0 II rsq == 0.0); 
fac = sqrt(-2.0*log(rsq)/rsq); 

/* Now make the Box-Muller transformation to get two normal deviates. */ 

/* Return one and save the other for the next call */ 

gset = vl*fac; 

iset =1; /* Set flag */ 

return v2*fac; 

} 

else 

{ 

/* We have an extra deviate handy. Unset the flag and return it */ 
iset = 0; 
return gset; 

} 


int bitgen(long *idum) 

{ 

double a; 
int bit; 

a = ran2(idum); 
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if (a >= 0.5) 
bit = 1; 
else 

bit = 0; 
return bit; 

} 


void add_noise(Iong *idum, double *data, double nn, double sigpow, double EbNo, 
double gain) 

{ 

int i; 

double var, sd, A; 

A = gain*pow(10.0, (EbNo/ 10.0)); 

var= sigpow*Nss/(2*A); 

sd = sqrt(var); 

for (i=0; icnn; i++) 
data[i] += sd*gasdev2(idum); 


B.3 Filter 

void filter(double *data, int M, double *h, double *x, long nn) 

{ 


int i, j; 

double sum, *x; 
x = dvector(2*M+l); 


for(i=0; i<=2*M; i++) x[i] = 0.0; 
x[2*M-l] = data[0]; 
x[2*M] = data[ 1 ] ; 

for(i=2; i<nn+2*M; i++) 

{ 

sum = 0.0; 
for(j~0; j<2*M; j++) 

{ 

sum += x[j]*h[j]; 

x{j] = xU+i]; 


x[2*M] = datafi]; 
data[i-2] = sum; 

} 

free(x); 
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B.4 Modulator/Demodulator 


void modulate (double *data, int *v , long N) 

{ 

int i, j; 
double A; 

A = sqrt(2.0*Eb/Tb); 
for(i=0;i<N;i++) 


if(v[i]=0) 


for(j=0;j<N$s;j++) 
data[i*Nss + j] = -A; 


else 

{ 


} 


for(j=0;j<Nss;j++) 
data[i*Nss + j] = A; 


void demod(doubie *data, int *out, long N) 

{ 

int i; 
double I; 

for(i=0;i<N;i++) 

{ 

1= data[i*Nss+Nss/2 -1] + data[i*Nss+Nss/2]; 

if(I >= 0.0) out[i] = 1; 
else out[i] = 0; 

} 


void demodsoft(double *data, int *out, long N, int **soft_metric, int Q, int a) 


{ 

int i, j, level, **num_bits_per_level; 
double I, b, amp, *prob; 

amp = 2.0; 

num_bits_per_level = int_matrix„2d(2, Q); 
prob = dvector(Q); 
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for(i=0; i<2; i++) 

{ 

for(j=0; j<Q; j++) 
num_bits_per_level[i][j] = 0; 

} 


for(i=0;i<N;i++) 

{ 

1= (data[i*Nss+Nss/2 -l] + data[i*Nss+Nss/2])/2.0; 
for(j=0; j<Q; j++) 

if(I>=amp*(2*j-Q)/Q && I<= amp*(2*(j+l)-Q)/Q) 
level = j; 

if( I >= amp) level = Q- 1 ; 
if(I <= -amp) level = 0; 

num_bits_per_level[out[i]][ level] += 1; 

out[i] = level; 


} 


for(i=0; i<Q; i++) 

prob[i] = (double)(num_bits_per_level[0][i] + num_bits_per_level[l][Q-i-l] + 1)/(2*N+Q); 

b = -log(prob[Q-l])/log(2.0); 

for(i— 0; i<Q; i++) 

{ 

soft_metric[0][i] = (int)(floor)(a*(log(prob[i])/log(2.0) + b)); 
soft_metric[l][Q-i-l] = soft_metric[0] [i] ; 

} 

free_2d_int_matrix(2, num_bits_per_level); 
free(prob); 


B.5 Calculate power 

double calc_power(double *data, long nn) 

{ 

int i; 

double sum, P; 

sum = 0.0; 
for(i=0; icnn; i++) 
sum += pow(data[i], 2); 

P = sum/(nn); 
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return P; 


} 


B.6 Interleaver and deinterleaver 

void interleave(int *v, int rows, int rs_n, int rs_m, int I, int num) 

{ 

/* rows = rs_m*rs_n , I = interleaving depth */ 

int ii, i, j, 1, **matrix; 

matrix = int_matrix_2d(rows, I); 

/* Fill matrix by column */ 

for (ii=0; iicnum; ii++) 

{ 

for(i=0; i< I; i++) 
for(j=0; j<rows; j++) 
matrix[j][i] = v[ii*rows*I + i*rows +jj; 

/* Exit matrix by row */ 

for (i=0; i<rs_n; i++) 
for(j=0; j<I; j++) 
for(l=0; l<rs_m; 1++) 

v[ii*rows*I + i*I*rs_m +j*rs_m + 1] = matrix[i*rs_m+l][j]; 


free_2d_int_matrix(rows, matrix); 


} 

void deinterleave(int *v, int rows, int rs_n, int rs_m, int I, int num) 

{ 

/* rows = rs_m*rs_n , I = interleaving depth */ 

int ii, i, j, 1, **matrix; 

matrix = int_matrix_2d(rows, I); 


/* Fill matrix by row */ 

for (ii=0; iicnum; ii++) 

{ 

for (i=0; i<rs„n; i++) 
for(j=0; jcl; j++) 
for(l=0; l<rs_m; 1++) 

matrix[i*rs_m+l][j] = v[ii*rows*I + i*I*rs_m +j*rs_m + 1] 
/* Exit matrix by column */ 
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for (i=0; i<I; i++) 
for(j=0; j<rows; j++) 
v[ii*rows*I + i*rows +j] = matrix[j][i]; 

} 


free_2d_int_matri x (ro ws , matri x) ; 


B.7 Convolutional encoder and decoders 

int bin2dec(int *temp, int a) 

{ 

int i, sum; 
sum = 0; 

for(i=0; i<a; i++) sum += temp[i]*pow(2, a-i-1); 
return sum; 

} 


void dec2bin( int *temp, int dec, int a) 

{ 

int i, sum, c; 

sum - dec; 
for(i=0; i<a; i++) 

{ 

c = pow(2, a-I-i); 
if (sum >= c) 

{ 

temp[i] = 1 ; 
sum -= c; 

} 

else temp[i] = 0; 

} 


void conv_encode( int *u, int n, int k, int m, int num) 

{ 

int ii, i, j, 1, K, truncjength, **g, **mem, *out; 
FILE *gen; 

gen = fopen("kl n2m6.dat", "r"); 

K = m+1; /* The constraint length */ 

truncjength = 100; 


g = int_matrix_2d(n, k*K); 
mem = int_matrix_2d(k, K); 
out= ivector(n*num+truncJength); 
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/* Obtain the generator matrix */ 


for(i=0; i<n; i++) 
for(j=0; j<K*k; j++) 
fscanf(gen, "%d", &g[i][j]); 

for(ii=0; iicnum+truncjength; ii++ ) 

{ 

/* Shift the contents of memory */ 
for(i=0; i<k; i++) 
for(j= K-2; j>=0; j— ) 
mem[i][j+l] = mem[i][j]; 

/* Insert new bits into the encoder *1 
for(i=0; i<k; i++) mem[i][0] = u[ii*k + i] ; 

/* Begin encoding process */ 
for(i=0; i<n; i++) 

{ 

out[ii*n + i] = 0; 
for(j=0; j<k; j++) 
for(l=0; 1<K; 1++) 

out[ii*n + i] g[i][j*K + 1] & mem[j][l]; 

} 


} 

for(i=0; i<n*num; i++) u[i] = out[i] ; 
free(out); 

free_2d_int_matrix(n, g); 
free_2d_int_matrix(k, mem); 
fclose(gen); 


int conv_decode(int *v, int n, int k, int m, int num, int **soft_metric) 


int iii, ii, i, j, 1, a, 11, K; 

int starting_state, partjnetric, old_state, flag; 

int max_metric, max_state, min_metric, num_states, num_inputs; 

int **g, **mem, **path, **path_next; 

int **prev_state, **branch_out, ^counter, **branchjn; 

int *metric, *prev_metric, *out, *state; 

int *input, !ic new_state, *r, truncjength; 

FILE *gen; 

truncjength = 32; 

K = m+1; /* The constraint length */ 
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num_states = pow(2, m*k); 
num_inputs = pow(2, k); 

prev_state = int_matrix_2d(num_states, num_inputs); 
branch_in = int_matrix_2d(num_states, num_inputs); 
branch_out = int_matrix_2d(num_states, numjnputs); 
path = int_matrix_2d(num_states, truncjength); 
path_next = int_matrix_2d(num_states, truncjength); 


g = int_matrix_2d(n, k*K); 
mem = int_matrix_2d(k, K); 
metric = ivector(num_states); 
prev_metric = ivector(num_states); 
counter = ivector(num_states); 
input = ivector(k); 
out = ivector(k); 
state = ivector(k*m); 
new_state = ivector(k*K); 
r = ivector(n); 

gen = fopen("kl n2m6.dat", "r"); 

/* Obtain the generator matrix */ 

for(i=0; i<n; i++) 
for(j=0; j<K*k; j++) 

{ 


fscanf(gen, "%d", &g[i][j]); 


fclose(gen); 

for(i=0; i<num__states; i++) counter[i] = 0; 

for(iii=0; iii< num_states; iii++) 

{ 

for(ii=0; ii< numjnputs; ii++) 

{ 

/* Obtain the binary representation of the state of the encoder */ 
dec2bin(state, iii, k*m); 

/* Initialize the encoder memory to the current state */ 
for(ll=0; ll<k; 11++) 
for(l-0; l<m; 1++) 
mem[ll][l] = state[il*m + 1]; 

/* Shift the contents of memory */ 
for(i=0; i<k; i++) 
for(j= K-2; j>=0; j— ) 
mem[i][j+l] = mem[i][j]; 
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/* Obtain the binary representation of the encoder input */ 
dec2bin(input, ii, k); 

/* Insert new bits into the encoder */ 
for(i-0; i<k; i++) mem[i][0] = input[i] ; 

/* Begin encoding process */ 
for(i=0; i<n; i++) 

{ 

out[i] = 0; 
for(j=0; j<k; j++) 
for(I=0; 1<K; 1++) 

out[i] A = g[i][j*K + 1] & mem[j][l]; 

} 


/* Find out the new state of the encoder */ 
for(ll=0; ll<k; 11++) 
for(l=0; l<m; 1++) 

new_state[ll*m+l] = mem[ll][l]; 


a = bin2dec(new_state, k*m); 

prev„state[a][counter[a]] - iii; 
branch_out[a] [counter [a]] = bin2dec(out,n); 
branch_in[iii][ii] = a; 

++counter[a]; 

} 

} /* End obtaining the decoder information */ 

starting_state = 0; 

for(i=0; i< num_states; i++) prev_metric[i] = -99; 
prev_metric[starting_state] = 0; 

for(iii=0; iii<num+ truncjength; iii++) 

{ 


/* Input the latest word */ 

for(i=0; i<n; i++) r[i] = v[iii*n+i]; 

for(i=0; i< num_states; i++) 

{ 


metric[i] = -99; 

for(j=0; j< num_inputs; j++) 

{ 

part_metric = prev_metric[prev_state[i][j]]; 
for(I=0; i<n; 1++) 

part_metric += soft_metric[((branch_out[i]jj] » l)&l)][r[n-l-l]]; 
if(part_metric > metric[i]) 
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{ 

metric[i] = part_metric; 
old_state = prev_state[i] [j] ; 

} 

} 

/* Now we have the old state and the metric, update the path info */ 

/* shift the path information */ 

for (1=0; 1< trunc_length - 1 ; 1++) 
path_next[i][l] = path[old_stale][l+l]; 

/* Insert the new branch input into the path */ 

for(l=0; l<num_inputs; 1++) 
if(branch_in[old_state][l] == i) 

{ 

path_next[i][trunc_length-l] = 1; 

} 

} 


/* Obtain the output for this iteration if trunc_length bits have entered the buffer */ 

if(iii >= truncjength- 1 ) 

{ 

max_metric = 0; 

/* Determine the path with the maximum metric */ 
flag = 0; 

for(i=0; i< num_states; i++) 

{ 

if(metric[i] > max_metric) 

{ 

max_metric = metricfi]; 
max_state = i; 

if(max_metric > 100000) flag = 1; 

} 

} 

if(flag) 

{ 

min_metric = max_metric; 
for(i=0; i< num_states; i++) 

{ 

if(metric[i] < min_metric) 

{ 

min_metric = metrical; 

} 

} 

for(i=0; i<num_states; i++) metric[i] -= min_metric; 

} 


138 



/* Obtain the output from the path with the maximum metric */ 


for(i=0; i<k; i++) v[(iii-trunc_length+2)*k-i-l] = ((path_next[max_state][0]»i)&l); 

} 


for(i=0; i< num_states; i++) prev_metric[i] = metricfi]; 

for(i=0; i<num_states; i++) 
for(j=0; j<trunc_length; j++) 
path[i][j] = path_next[i][j]; 


free_2d_int_matrix(num_states, prev_state); 

free_2d_int_matrix(num_states, branchjn); 

free_2d_int_matrix(num_states, branch_out); 

free_2d_int_matrix(num_states, path); 

free_2d_int_matrix(num_states, path_next); 

free_2d_int_matrix(n, g); 

free_2d_int_matrix(k, mem); 

free(metric); 

free(prev_metric); 

free (counter); 

free(input); 

free(state); 

free(new_state); 

free(r); 

free(out); 

return(O); 


} 


int SOVA(int *v, int n, int k, int m, int num, int **soft_metric, int al) 

{ 


int iii, ii, i, j, 1, a, II, K, s 1 , s2, m2; 

int starting_state, part_metric[2], old_state[2], flag; 

long max_metric, max_state, min_metric, num_states, num_inputs; 

int **g, **mem, **path, Ht *path_next,**prev_state; 

int **branch_out, *conv_counter, **branch_in, *metric; 

int *prev_metric, *out, *state, *input, *new_state, *r; 

double c, delta, **L, **L_next; 

int L_Q, L_max, nlevels, truncjength; 
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FILE *gen; 


gen = fopen("kl n2m6.dat", "r"); 

K = m+1; /* The constraint length */ 

L_Q = 8; 

L_max = 8; 
truncjength = 32; 
num_states = pow(2, m*k); 
num_inputs = pow(2,k); 
nlevels = pow(2, L_Q); 

g = int_matrix_2d(n, k*K); 
mem = int_matrix_2d(k, K); 
metric = ivector(num_states); 
prev_metric = ivector(num_states); 
conv_counter = ivector(num_states); 
input = ivector(k); 
out = ivector(n); 
state = ivector(k*m); 
new_state = ivector(k*K); 
r = ivector(n); 

prev_state = int_matrix_2d(num_states, num_inputs); 
branch_in = int_matrix_2d(num_states, num_inputs); 
branch_out = int_matrix_2d(num_states, num_inputs); 
path = int_matrix_2d(num_states, trunc_length); 
path_next = int_matrix_2d(num_states, truncjength); 
L = double_matrix_2d(num_states, truncjength); 
L_next = double_matrix_2d(num_states, trunc_length); 


/* Obtain the generator matrix */ 

for(i=0; i<n; i++) 
for(j=0; j<K*k; j++) 

{ 

fscanf(gen, "%d", &g[i][j]); 

} 

fclose(gen); 

num_states = pow(2, m*k); 
num_inputs - pow(2, k); 

c = log(2.0)/(double)(al); 

for(i=0; i<num„states; i++) conv_counter[i] = 0; 

for(iii=0; iii< num_states; iii++) 

{ 

for(ii=0; ii< num_inputs; ii++) 

{ 
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/* Obtain the binary representation of the state of the encoder */ 
dec2bin(state, iii, k*m); 


/* Initialize the encoder memory to the current state */ 
for(ll=0; lkk; 11++) 
for(l=0; l<m; 1++) 
mem[ll][l] = state[ll*m + 1]; 


/* Shift the contents of memory */ 
for(i=0; i<k; i++) 
for(j- K-2; j>=0; j — ) 
mem[i][j+l] = mem[i][j]; 


/* Obtain the binary representation of the encoder input */ 
dec2bin(input, ii, k); 

/* Insert new bits into the encoder */ 
for(i=0; i<k; i++) mem[i][0] = input[i]; 

/* Begin encoding process */ 
for(i=0; i<n; i++) 

{ 

out[i] = 0; 
for(i=0; j<k; j++) 
for(l=0; 1<K; 1++) 

out[i] A = g[i][j*K + 1] & mem[j][l]; 

} 


/* Find out the new state of the encoder */ 
for(ll=0; lkk; 11++) 
for(l=0; km; 1++) 

new_state[ll*m+l] = mem[ll][l]; 


a - bin2dec(new_state, k*m); 

prev_state[a][conv_counter[a]] = iii; 
branch_out[a][conv_counter[a]] = bin2dec(out,n); 
branch_in[iii][ii] = a; 

++conv_counter[a] ; 

} 

} /* End obtaining the decoder information */ 

starting_state = 0; 

for(i=0; k num_states; i++) prev_metric[i] = -99; 
prev_metric[starting_state] = 0; 

for(i=0; i<num_states; i++) 
for(j=0;j<trunc_length; j++) 

L[i][j] = 99999.0; 


141 



for(iii=0; iii<num+trunc_Jength; iii++) 

{ 


/* Input the latest word */ 


for(i=0; i<n; i++) r[i] = v[iii*n+i]; 


for(i=0; i< num_states; i++) 

{ 


metric[i] = -99; 
m2 = -99; 

for(j=0; j< num_inputs; j++) 

{ 

old_state[j] = prev__state[i][j] ; 
part_metric[j] = prev_metric[old_state[j]]; 
for(l=0; l<n; 1++) 

part_metric[j] += soft_metric[((branch_out[i][j] » l)&l)][r[n-l-l]]; 


/* Determine the maximum metric and the second maximum metric */ 

for(j=0; j<num_inputs; j++) 

{ 

if(part„metric[j] >= metric[i]) 

{ 

si = old_state|j]; 
metric[i] = part_metric[j]; 

a = j; 

} 

} 

for(j=0; j<num_inputs; j++) 

{ 

if G==a) continue; 
if(part_metric[j] >= m2) 

{ 

s2 = old_state[j]; 
m2 = part_metric[j]; 

} 


/* Now we have the old state and the metric, update the path info */ 

/* Calculate delta */ 

delta = (double)c*(metric[i] - m2); 

/* Update the reliability information */ 

for(j=l; j<=trunc_length-m+l ; j++) 

{ 

if(path[sl][j] != path[s2][j]) 
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{ 

if(L[s 1 ][j] >= delta) 
L_next[i][j-1] = delta; 
else 

L_next[i][j-1] =L[sl][j]; 


else 

L_next[i] [j- 1 ] - L[s 1 ] [j] ; 


for(j=trunc_length-m+2; j<trunc_length; j++) 
L_next[i][j-1] = L[s 1 ][j] ; 


/* shift the path information */ 

for (1=0; l<trunc_length- 1 ; 1++) 
path_next[i][l] = path[s !][!+!]; 


/* Insert the new branch input into the path */ 

for(l=0; l<num_inputs; 1++) 
if(branch_in[sl][l] == i) 

{ 

path„next[i][trunc_length-l] = 1; 

} 

L„next[i][trunc_length-1] = 10000.0; 


/* Obtain the output for this iteration if truncjength bits have entered the buffer */ 

if(L_DEBUG) 

{ 

printf("\nReliabilities\n"); 

for(l=0; l<num_states; 1++) 

{ 

printf("\n”); 

for(i=0; i<trunc_length; i++) 

{ 

printf("%3.31f ", L_next[I][i]); 

} 

} 

if(PAUSE) getchar(); 

printf("\nbits in paths\n"); 

for(l=0; l<num_states; 1++) 

{ 

printf("\n"); 
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for(i=0; i<trunc_length; i++) 

{ 

printf("%d ", path_next[l][i]); 

} 

} 

if(PAUSE) getchar(); 

} 

if(iii >- trunc_length-l) 

{ 

if((iii-trunc_length+l)%8 == 0) 

{ 

max_metric = 0; 

/* Determine the path with the maximum metric */ 
flag = 0; 

for(i=0; i< num_states; i++) 

{ 

if(metric[i] > max_metric) 

{ 

max_metric = metric[i]; 
max_state = i; 

if(max_metric > 1000000) flag = 1; 

} 

} 


if(flag) 

{ 

min_metric = max_metric; 
for(i=0; i< num_states; i++) 

{ 

if(metric[i] < min_metric) 

{ 

min_metric = metric[i]; 

} 

} 

for(i=0; i<num_states; i++) metric[i] -= min_metric; 

} 

/* Obtain the output from the path with the maximum metric */ 

for(j=0; j<8; j++) 

{ 

for(i=l; i<= nlevels; i++) 

{ 

if((L_next[max_state][j] >= (double)(i*L_max/nlevels)) && 
(L_next[max_state][j] < (double)((i+I)*L_max)/nlevels)) 

{ 

if(path_next[max_state][j] == 1) v[iii-trunc_length+l+j] = i 
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else v[iii-trunc_length+ l+j] = -i; 


break; 


} 

} 

if(L_next[max_state][j] > L_max) 


if(path_next[max_state][j] == 1) v[iii-truncjength+l+j] = nlevels; 
else v[iii-truncjength+l+j] = -nlevels; 

} 


if(L_next[max_stateJ01 < L_max/nlevels) 

{ 

if(path_next[max_state][j] == 1) v[iii-trunc_length+l+j] = 1; 
else vfiii-truncjength+l-f-j] = -1; 



for(i=0; i< num_states; i++) prev_metric[i] = metric [i] ; 

for(i=0; i<num_states; i++) 
for(j=0; j<trunc„length; j++) 

{ 

path[i][j] = path_next[i][j]; 

L[i][j] = L_next[i][j]; 

} 


free(metric); 

free(prev_metric); 

free(conv_counter); 

free(input); 

free(out); 

free(state); 

free(new_state); 

free(r); 

free__2d_int_matrix(n, g); 
free_2d_int_matrix(k, mem); 
free_2d_int_matrix(num_states, prev_state); 
free_2d_int_matrix(num_states, branch_in); 
free_2d_int_matrix(num_states, branch_out); 
free„2d_int_matrix(num_states, path); 
free_2d_int_matrix(num_states, path_next); 
free_2d_double_matrix(num_states, L); 
free_2d_doubIe_matrix(num_states, L_next); 

return(O); 

} 
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B.8 Reed Solomon encoders and decoders 

void make_GF_table(int m, int *GF_table, int *dec_table, int GF_poly) 

/* This procedure generates the Galois Field GF(2 A m). Information */ 
/* is stored in two lookup tables (decjable and GF_table). */ 

/* */ 
/* GF_table[i] = base ten number, where 'i' is the power of alpha */ 
/* dec_table[i] = power of alpha, where 'i' is the base ten repersentation. */ 
/* 

/* For Example: InGF(16) */ 

/* */ 

/* a A 6 = 0 0 11 = 3 (Base ten, decimal equivalent) */ 

/* */ 

/* is the same as */ 

/* */ 

/* GF_talble[6] - 3 or dec_table[3] = 6 */ 

/* */ 

/* These two tables are used to transform a number back and fourth */ 

/* between the two representations (GF and decimal). */ 

/* This is done in order to perform Galois Field arithmatic. When */ 

/* adding, it is easier to use the decimal representation, because adding */ 
/* two numbers, a and b, is a A b, where ' A ' is exclusive-or. */ 

/* If a and b were in GF powers of alpha representation, then adding */ 
/* the two numbers becomes: */ 

/* */ 

/* GF_table[a] A GF_table[b] = c */ 

/* */ 

/* Multiplying two numbers is easier to use the powers of alpha */ 

/* representation. For example, inGF(16) */ 

/* alpha A 4 x alpha A 5 = alpha A 9 */ 

/* Multiplication is acheived by simply adding the powers of alpha (4+5)*/ 
/* This addition is modulus q-1. */ 

/* alpha A 14 x alpha A 7 = alpha A 21 mod (15) = alpha A 7 */ 

/* Due to the fact that the element zero is not represented in the */ 

/* GF_table because in the powers of alpha representation, 0 is actually */ 
/* equal to alpha A 0 = 1. So the GF representation for the element 0 is -1 */ 
/* Care has to be taken when multiplying two elements. The multiplying */ 
/* can only take place if a or b is not equal to - 1 . Else the result is equal */ 
/* to zero. So dec_table[0] = - 1 , but GF_table[-l] does not equal */ 
/* zero. This returns an error because indixes of arrays can not have */ 
/* negative values. */ 

/* */ 

/* GF_poly is the irriducable (or primitave) polynomial that is used to */ 
/* generate the field. For example, for m = 4, the irriducable */ 

/* polynomial is: */ 

/* */ 

/* X A 4 + X + 1 = 0 or X A 4 = X+ 1 = 1100 = 12 */ 

/* */ 

/* The X A 4 = 12 is the way GF_poly is stored. The field is created by */ 
/* repeated right shifts (multiplication by alpha). If a T is right shifted */ 
/* from the right most position (corresponding to alpha A 3), we will have */ 
/* alpha A 4. But alpha A 4 is equal to alpha +1(1 100 = 12), so 12 is */ 
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/* 'x-or' ed with the result. 


*/ 


{ 

int i, q; 

GF_table[0] = pow(2, (m-1)); 
q = pow(2, m); 


for( i=l; i<q-l; i++) 

{ 

if (GF_table[i- 1 ] & 1 =- 1) 

{ 

GF_table[i] = GF_table[i-l] » 1; 
GF_table[i] = GF_tabIe[i] A GF_poly; 

} 

else GF_table[i] = GF_table[i- 1 ] » 1; 

} 

for(i=0; i<q-l; i++) dec_table[GF_table[i]] = i; 
dec_table[0] = - 1 ; 

} 


void rs_encode (int *data, int n, int k, int m, int num) 


/* The following procedure encodes a data vector using a Reed Solomon */ 
/* code. The parameters are defined as follows: */ 

/* */ 
/* m = number of data bits per symbol */ 

I* q = 2 A m specifies the field GF(q) = GF(2 A m) */ 

/* n = q- 1 the block length */ 

/* t = maximum number of errors that can be corrected */ 

/* k = n - 2*t the number of information bits per block */ 

/* num = the total number of blocks to be encoded */ 

/* data = the uncoded data vector */ 

/* */ 
/* The field is first generated by calling the procedure 'make_GF_table\ */ 
/* Then the systematic encoding is accomplished by using feedback shift */ 
/* registers to generate the parity symbols. The connections are specified */ 
/* by the generator polynomial g[]. */ 


int GF_poly[l 1] = {0,0,0,6,12,20,48,72,184,272,576}; 
int i, ii, j, 1, cl, t, sum; 

int *u, *vv, *g, *datal, *GF_table, *dec_table; 
t= (n-k)/2; 

GF_table = ivector(n+l); 

dec_table = ivector(n+l); 

g = ivector(n-k+l); 

u = ivector(k); 

vv = ivector(n); 

datal = ivector(n*m*num); 
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make_GF_table (m, GF_table, dec_table, GF_poly[m]); 


for (i=0; i<= n-k; i++) g[i] = 0; 


g[0] = GF_table[l]; 
g[l] = GF_table[0]; 

/* Create the generator polynomial g(X) */ 
for (i=2; i<= n-k; i++) 

{ 

g[i] = GF_table[0]; 
for (j— i- 1 ; j>0; j — ) 

if (g[j] != 0) g[j] = g[j-l] A GF_table[ (dec_table[g[j]] + i) % n]; 
else g|j] = glj-1]; 

g[0] = GF_table[ (dec_table[g[0]] + i) % n] ; 

} 


/* Change g to GF representation */ 
for(i=0; i<= n-k; i++) g(i] = dec_table[g[i]]; 

for (ii = 0; ii <num; ii++) 

{ 

for(j = 0;j<k; j++) 

{ 

u[j] = 0; 

for (1 = 0; 1 < m; 1++) 

u[j] += data[ii*k*m + j*m +l]*pow(2, m-1-1); 

1 

for (i=0; i< n-k; i++) vv[i] = 0; 
for (i=k- 1 ; i>=0; i— ) 


cl = dec_table[u[i] A vv[2*t-l]]; 
if (cl != -1) 

{ 

for(j=2*t- 1 ; j>0; j— ) 

if (g[j] ■= -0 VV U] = vvfj- 1 ] A GF_table[(g[j]+c 1 )%n]; 
else vv[j] = vv[j-l]; 
vv[0] = GF_table[(g[0] + cl)%n]; 

} 

else 


for (j= n-k- 1 ; j>0; j — ) 
vv[j] = vv[j-l]; 
vv[0] = 0; 

} 

} 

for(i=0; i<k; i++) vv[i+n-k] = u[i] ; 

for(i=0; i<n; i++) 

{ 

sum = vv[i]; 
for (j=0; j<m; j++) 

{ 

cl = pow(2, m-j-1); 
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if (sum >- cl) 

{ 

datal [ii*n*m + i*m + j] = !; 
sum -= c l ; 

} 

else datal [ii*n*m + i*m + j] = 0; 

} 

} 

} 

for (i=0; i<n*m*num; i++) datafi] = datal [i]; 

free(datal); 

free(GF_table); 

free(dec_table); 

free(vv); 

free(g); 

free(u); 

} 

void rs_decode(int *data, int ^counter, int n, int k, int q, int m, int num) 

/* The following procedure decodes a Reed Solomon encoded data */ 

/* vector. This RS decoder has been modified to handle erasures. An */ 
/* erased position is denoted '-2'. The code can correct e errors and f */ 
/* erasures if 2*e + f <= dmin where dmin is equal to n - k or 2*t. */ 

/* The parameters are defined as follows: */ 

/* */ 

/* m = number of data bits per symbol */ 

/* q = 2 A m specifies the field GF(q) - GF(2 A m) */ 

/* n = q-l the block length */ 

/* t = maximum number of errors that can be corrected */ 

/* k = n - 2*t the number of information bits per block */ 

/* num = the total number of blocks to be encoded */ 

/* data = the uncoded data vector */ 

/* lambda = The error locator polynomial */ 

/* S = The syndrome polynomial */ 

/* */ 

/* The decoding is done using the Berlekamp-Massey algorithm. Details */ 
/* of the algorithm can be found in "Theory and Practice of Error Control */ 
/* Codes" by Richard E. Blahut. This procedure can handle non-erasure */ 
/* decoding also. */ 

{ 

int GF_poly[ll]= {0,0,0,6,12,20,48,72,184,272,576}; 

int i, ii, j, 1, t, rr, sum, L, cl , degjambda, num_erasure, delta_r; 

int erasure, error, decode_fiag, num_errors, numer, den; 

int *GF_table, *dec_table, *r, *U, *S, *lambda, *omega; 

int *B, *T, *tmp, *beta; 

t = (n-k)/2; 
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GF_table = ivector(q); 
dec_table = ivector(q); 
beta = ivector(n); 
r = ivector(n); 

U = ivector(n); 

S - ivector(n); 

T = ivector(n); 

B = ivector(n); 
tmp = ivector(n); 
lambda = ivector(n); 
omega = ivector(2*t+l); 


/* Create GF(2 A m) field */ 

make_GF_table (m, GF_table, dec_table, GF_poly[m]); 


for (ii = 0; ii < num; ii ++) 

{ 


decode_flag = 0; 
num_erasure = 0; 
for(j = 0; j < n; j++) 

{ 

r[j] = 0; 

sum = 0; 

erasure = 0; 

for (1 = 0; 1 < m; 1++) 

{ 

if (data[ii*n*m + j*m +1] == -2) /* Erased bit */ 

{ 

U[num_erasure] = j; 
if (erasure — 0) 

{ 

num_erasure++; 
erasure = 1 ; 

} 

rUl = -i; 

datafii^n^m + j*m +1] = 0; 
break; 

} 

else sum += data[ii*n*m + j *m +1] *pow(2, m- 1-1); 

} 

if (r[j] != -1) r[j] = dec_table[sum]; 

} 

error = 0; 

/* Compute the syndrome *1 
for(i=l; i<= n-k; i++) 


S[i] = 0; 

for(j=0; j<n; j++) 
if(r[j] !=-D 

S[i] A = GF_table[ (r[j] + i*j) % n]; 
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if (S[i] != 0) error = 1 ; /* If nonzero syndrome, there is an error*/ 

} 

S[0] =0; 

/* Convert to GF representation */ 
for(i=l ; i<=n-k; i++) S[i] = decjable[S[i]]; 


if (error) /* If the syndrome is equal to zero, no decoding nessasary */ 

{ 

for(i=0; i < n; i++) 

{ 

lambda[i] = 0; 

B(i] = -1; 

T[i] = 0; 

} 

lambda[0] = GF_table[0]; /* =1 */ 

L = 0; 

deg_lambda = 0; 

B[0] = 0; /* = 1 */ 

for(rr=l; rr <= 2*t; rr++) 

{ 

if (rr <= num_erasure) 

{ 

for (i=l; i <= deg Jambda+ 1 ; i++) 

{ 

if (lambda[i-l] != 0) 

tmp[i] = GF_table[ (U[rr-1] + dec_table[ lambda[i-lj ] ) % n]; 
else tmp[i] = 0; 

} 

for (j== 1 ; j <= deg Jambda+ 1 ; j++) lambda[j] A = tmp[j]; 
deg_lambda++; 

for(j=0; j<= degjambda; j++) B [j] = dec_table[lambda[j]] ; 

L ++j 

} 


{ 

/* Compute the discrepancy */ 
delta_r = 0; 
for(j=0; j<=rr; j++) 
if(lambda[j] != 0 && S[rr-j] !=-l) 
delta_r A = GF_table[ (dec_table[lambda[j]] + S[rr-j]) % n]; 
delta_r = dec_table[delta_r]; 

if (delta_r != -1) 

{ 

/* T(x) < — lambda(x) + delta_r * x * B(x) */ 
for (j=l ; j <= degjambda + 1 ; j++) 

if (B[j- 1 ) !- -1) 

T[j] = lambda[j] A GF_table[ (delta_r + B[j- 1 ]) % n]; 

T[0] = lambda[0]; 
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++deg_lambda ; 


if (2*L <= rr + num_erasure -1) 

{ 

L = rr - L + num_erasure; 

/* B(x) <— B(x) / delta_r */ 
for(i=0; i<= degjambda; j++) 
if( lambdafj] != 0) 

B[j] = (n-delta_r + dec_table[lambda[j]]) % n ; 
else B[j] - - 1 ; 

/* lambda(x) < — T(x) */ 

for(j=0; j<=deg_lambda+ 1 ; j++) lambda[j] = T[j]; 


} 

else 

I 

/* lambda(x) < — T(x) */ 
for(j=0; j<n-k; j++) lambda[j] = T[j] ; 

/* B(x) < x * B(x) */ 

tmp[0] = -1; 

for(j=l; j<=n-k; j++) tmpfi] = B [j- 1 ] ; 
for(j=0; j<=n-k; j++) B[j] = tmpU]; 

} 


else 

{ 

/* B(x) < x * B(x) */ 

tmp[0] = -1; 

for(j= 1 ; j<=n-k; j++) tmp[j] = B[j-1]; 
for(j=0; j<=n-k; j++) B[j] = tmp[j]; 

} 

} 

} 


/* Change lambda(x) to GF representation */ 
for(i=0; i<n; i++) lambda[i] = dec„table[lambda[i]]; 

/* Compute the degree of lambda(x) */ 
deg_lambda = n; 
for(i=n-l; i>=0; i— ) 

if (lambda[i] != -1 && deg_lambda == n) deg_lambda = i; 

if (deg_lambda <= 2*t) /* Below the capacity of the code */ 

{ 


/* Comupte omega(X) = [l+S(X)j * lambda(X) */ 

for(i=0; i<=n-k; i++) omega[i] = 0; 
for (i=0; i<=n-k; i++) 


for (j=0; j<=n-k; j++) 
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if ((i+j) >= n-k+1) continue; 

if (S[i] != -1 && lambda[j] != -1) 

omega[i+j] A = GF_table[ (S[i] + lambda[j]) % n ]; 


/* Convert omega(x) to GF representation */ 

for (i=0; i<=n-k; i++) omega[i] = dec_table[omega[i]]; 

/* Find the roots of Iambda(X). The inverses of the roots gives us the */ 
/* location of the errors. */ 

num_errors = 0; 
for (i=0; i<q-l; i++) 

{ 

sum = 0; 

for (j=0; j<2*t; j++) 
if (lambdafj] != -1) 

sum A = GF_table[ (lambda[j] + i*j) % (q-1) ]; 
if (sum == 0) 

{ 

beta[num_errors] = (n-i)%n; 
num_errors++; 

} 

} 


if( (2*(num_errors-num_erasure) + num_erasure) <= 2*t) 

{ 

if (num_errors == degjambda) 

{ 


/* Convert r to base 10 representation */ 

for(i=0; i<n; i++) 
if (r[i] != -1) r[i] = GF_table[r[i]]; 
else r[i] = 0; 


/* Calculate the error values and correct the received vector */ 

for (i=0; i<num_errors; i++) 

{ 

/* Calculate the denominator */ 
den = 0; 

for (j=0; j<=2*t; j++) 

{ 

if (j%2 == 0) continue; 
if (lambdaQ] != -1) 

{ 

cl = GF_table[ ( (q-l-beta[i])*(j-l) + lambdafj]) % n ]; 
den A = c 1 ; 

} 

} 
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den = dec_table[den]; 

/* Calculate the numerator */ 

numer = 0; 

for (j=0; j<=2*t; j++) 

if (omega[j] != -1) numer A = GF_table[ ((q-l-beta[i])*j + omega[j])%n ]; 

numer = dec_table[numer]; 

/* Correct the erroneous value */ 
if(numer != -1) 

r[beta[i]] A = GF_table[ (n + numer + beta[i] - den) % n]; 


/* Change r back into the GF representation */ 

for(i=0; i<n; i++) r[i] = dec_table[r[i]]; 

} /* end if (num_errors == deg_lambda) */ 

else 

{ 

decode_flag = 1 ; 
counter[3] += 1 ; 

/*printf("\n3"); 

*/ 

} 

} /* end if (num_erasure + num errors <= dmin) */ 

else 

{ 

decode_flag = 1 ; 
counter[2] += 1 ; 

/*printf("\n2"); 

*/ 

} 

} /* end if (deg_lambda <= 2*t) */ 

else 

{ 

decode_flag - 1 ; 
counter! 1] += 1; 

/*printf("\nl"); 

*/ 

) 

} /* end if (error) */ 

else 

{ 

decode_flag = 0; 
counter[0] += 1; 

/*printf("\nl"); 

*/ 

} 

if (decode^flag) 
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{ 

counter[5] += 1; 

/*printf("\n5"); 

*/ 

for(j=0; j<k; j++) 

{ 

for(l=0; l<m; 1++) 

data[ii*k*m + j*m + 1] = data[ii*n*m + (n-k+j)*m + 1]; 

} 


else 

{ 

counter[4] += 1; 

/*printf("\n4"); 

*/ 

for(j=0 ; j<k; j++) 

{ 

if (r[n-k+j] == -1) sum = 0; 
else sum = GF_table[r[n-k+j]l; 
for (1=0; l<m; 1++) 

{ 

cl = pow(2, m-1-1); 
if (sum >= cl) 

{ 

data [ii*k*m + j*m + 1] = 1; 
sum -= c 1 ; 

} 

else data [ii*k*m + j*m + 1] = 0; 



free(GFjable); 

free(dec_table); 

free(beta); 

free(r); 

ffee(U); 

free(S); 

free(T); 

free(B); 

free(tmp); 

free(lambda); 

free(omega); 


void rs_decode_erasure_method_l(int *data, int n, int k, int q, int m, int num, 
int *erasure_counter, int ^counter) 

/* The Reed Solomon deocder is modified to use erasure Method 1. This is to be */ 
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/* used with the SOVA. The bit reliabilities are first converted into symbol */ 

/* reliabilities. A table of the least reliabile symbols are compiled for each */ 

/* codeword. Decoding is attempted using no erasures. If successful, decoding */ 
/* stops. If unsuccessful, then decoding is reattemped using 2 erasures. For each */ 
/* unsuccessful decoding trial, two more erasures are added until either a */ 

/* successful decoding takes place, or a maximum number of erasures has been */ 
/* reached (16 for this simulation) */ 


{ 

int GF_poly[ll] = {0,0,0,6,12,20,48,72,184,272,576}; 
int i, ii, j, 1, t, rr, sum, L, cl, deg_lambda, delta_r; 
int error, decode_flag, num_errors, numer, den; 
int num_erasures, temp, max_eras, Reliability; 
int *GF_table, *dec_table, *r, *U, *S; 
int *lambda, *omega, *B, *T, *tmp, *beta, *min; 

t = (n-k)/2; 

max_eras = t; 


GF_table = ivector(n+l); 
dec_table = ivector(n+l); 
beta = ivector(n); 
r = ivector(n); 

U = ivector(2*t); 

S = ivector(n); 

T = ivector(n); 

B = ivector(n); 
tmp = ivector(n); 
lambda = ivector(n); 
omega = ivector(2*t+l); 
reliability = ivector(n); 
min = ivector(2*t); 


/* Create GF(2 A m) field */ 

make_GF_table (m, GF_table, dec_table, GF_poly[m]); 


for (ii = 0; ii < num; ii ++) 

{ 

for(j = 0; j < n; j++) 

{ 

r[j] = 0; 

reliability [j] = 500; 
sum = 0; 

for (1 = 0; 1 < m; I++) 

{ 

if( reliability [j] > abs(data[ii*n*m + j*m +1])) 
reliability [j] = abs(data[ii*n*m + j*m +1]); 

if(data[ii*n*m + j*m +1] <= 0) 
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data[ii*n*m + j*m +1] = 0; 
else 

data[ii*n*m + j*m +1] = 1; 
sum += data[ii*n*m + j*m +l]*pow(2, m-1-1); 

} 

r[j] = dec_tabie[sum]; 


/* Find the minimum reliabilities */ 

for(i=0; i<max_eras; i++) 

{ 

min[i] = reliability[i] ; 

U[i] = i; 

} 

/* Do the initial sorting. Place the minimum in min[0] and the */ 
/* maximum in min [max_eras - 1] */ 

for(i=0; i<max_eras; i++) 

{ 

for(j=0; j<max_eras - 1 ; j++) 

{ 

if(min[j] > min[j+l]) 

{ 

temp = min[j]; 
min[j] = minfj+1]; 
min[j+l] = temp; 
temp = U[j]; 

U[j] = U[j+I ]; 

UU+1] = temp; 

} 


} 


for(i=max_eras; i<n; i++) 

{ 

if(reliability[i] < min[max_eras - 1]) 

{ 

min[max_eras-l] = reliability [i]; 
U[max_eras-1] = i; 
for(j=max_eras- 1 ; j>= 1 ; j-) 

{ 

if(min[j] < min[j-l]) 

{ 

temp = min[j]; 
min[j] = min[j-l]; 
min[j-l] = temp; 
temp = U[j]; 

U[j] = U[j-13; 

U[j-13 = temp; 

} 
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else break; 


} 


} 

} 

/* The symbols with the minimum reliabilities are in positions */ 

/* U[0] ... U[num_erasures-1]. These symbols will be erased. */ 
/* Set these symbols equal to 0 before computing the syndrome */ 


decode_flag = 1 ; 

for(num_erasures=0; num_erasures <= max_eras; num_erasures +=2) 

{ 

if(decode_flag) 

{ 

for(i=0; i<num_erasures; i++) r[U[i]] = -l; /* =0 */ 
error = 0; 

/* Compute the syndrome */ 
for(i=l; i<= n-k; i++) 

{ 

S[i] = 0; 

for(j=0;j<n; j++) 

if(r[j] !=-l) 

S[i] A = GF_table[ (r[j] + i*j) % n]; 
if (S[i] != 0) error = 1 ; /* If nonzero syndrome, there is an error*/ 

} 

S[0] = 0; 

/* Convert to GF representation */ 
for(i=I; i<=n-k; i++) S[i] = dec_tabIe[S[i]]; 

if (error) /* If the syndrome is equal to zero, no decoding nessasary */ 

{ 

for(i=0; i < n; i++) 

{ 

lambda[i] = 0; 

B[i] = -1; 

T[i] = 0; 

} 

lambda[0] = GF_table[0]; /*=[*/ 

L = 0; 

deg_lambda = 0; 

B[0] = 0; /* = 1 */ 

for(rr=l; rr <= 2*t; rr++) 

{ 

if (rr <= num_erasures) 
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{ 

for (i= 1 ; i <= deg_lambda+ 1 ; i++) 

{ 

if (Iambda[i-i] != 0) 

tmp[i] = GF_table[(U[rr-l]+dec_table[lambda[i-l]])%n]; 
else tmp[i] = 0; 

} 

for (j^ 1 ; j <= degjambda+l; j++) lambda[j] A = tmp[j]; 
degjambda++; 

for(j=0; j<= degjambda; j++) B[j] = dec_table[lambda[j]]; 
L ++; 

} 

else 

{ 

/* Compute the discrepancy */ 
delta_r = 0; 
for(j=0; j<=rr; j++) 
if (lambda[j] !=0 && S[rr-j] !=-l) 
delta_r A = GF_table[(dec_table[lambda[j]]+S[rr-j])%n]; 
delta_r = dec_table[delta_r]; 

if (deltajr != -1) 

{ 

/* T(x) <— lambda(x) + delta_r * x * B(x) */ 

for (j=l; j <= degjambda +1; j++) 
if (B[j- 1 ] != -1) 

T[j] = lambdajj] A GFjable[ (delta_r + B[j- 1 ]) % n]; 
T[0] = lambda[0]; 

++deg_lambda; 

if (2*L <= rr + num_erasures-l) 

{ 

L = rr - L + num_erasures; 

/* B(x) < — B(x) / delta_r */ 
for(j=0; j<= degjambda; j++) 
if( lambda[j] != 0) 

B[j] = (n-delta_r+decjable[lambda[j]])%n; 
else B [j] = - 1 ; 

/* lambda(x) < — T(x) */ 

for(j=0; j<=degjambda+ 1 ; j++) lambda[j] = T[jj; 


} 

else 

{ 

/* lambda(x) <- — T(x) */ 
for(j=0; j<n-k; j++) lambda[j] = T[j] ; 

/* B(x) < x * B(x) */ 

tmp[0] = -1; 

for(j= 1 ; j<=n-k; j++) tmp[j] = B[j-1]; 
for(j=0; j<=n-k; j++) BJj] = tmp[j]; 


159 



} 


else 

{ 

/* B(x) < x * B(x) */ 

tmp[0] = -1; 

for(j=l;j<=n-k; j++) tmp[j] = B[j-1]; 
for(j=0; j<=n-k; j++) B[j] = tmp[j]; 

} 

} 

} 


/* Change lambda(x) to GF representation */ 
for(i=0; i<n; i++) lambda[i] = dec_table[lambda[i]]; 

/* Compute the degree of lambda(x) */ 
deg_lambda = n; 
for(i=n-l; i>=0; i--) 

if (lambda[i] != -1 && deg_lambda == n) deg_lambda = i; 

if (deg_!ambda <= 2*t) /* Below the capacity of the code */ 

{ 

/* Comupte omega(X) = [l+S(X)j * lambda(X) */ 

for(i=0; i<=n-k; i++) omega[i] = 0; 
for (i=0; i<=n-k; i++) 


for (j=0; j<=n-k; j++) 

{ 

if ((i+j) >= n-k+1) continue; 

if (S[i] != -1 && lambdafj] != -1) 

omega[i+j] A = GF_table[ (S[i] + lambda[j]) % n ]; 

} 

} 


/* Convert omega(x) to GF representation */ 

for (i=0; i<=n-k; i++) omega[i] = dec_table[omega[i]]; 

/* Find the roots of lambda(X). The inverses of the roots gives us the */ 
/* location of the errors. */ 

num_errors = 0; 
for (i=0; i<q- 1 ; i++) 

{ 

sum = 0; 

for (j=0; j<2*t; j++) 
if (lambdajj] != -1) 

sum A = GF_table[ (lambdaQ] + i*j) % (q-1) ]; 
if (sum == 0) 

{ 

beta[num_errors] = (n-i)%n; 
num_errors++; 

} 
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} 

if( (2*(num_errors-num_erasures) + num_erasures) <= 2*t) 

{ 

if (num_errors == degjambda) 

{ 


/* Correct Decoding, record the umber of erasures required */ 
erasure„counter[num_erasures] += 1 ; 
decode_flag = 0; 

/* Convert r to base 10 representation */ 

for(i=0; i<n; i++) 
if (r[i] !- -1) r[i] = GF_table[r[i]]; 
else r[i] = 0; 

/* Calculate the error values and correct the received vector */ 

for (i=0; i<num_errors; i++) 

{ 

/* Calculate the denominator */ 
den = 0; 

for (j=0; j<=2*t; j++) 

{ 

if (j%2 == 0) continue; 
if (lambdalj] != -1) 

{ 

cl = GF_table[ ( (q- 1 -beta[i])*(j- 1 ) + lambda[j]) % n ]; 
den A = c 1 ; 

} 

} 

den = dec_table[den]; 

/* Calculate the numerator */ 

numer = 0; 

for (j=0; j<=2*t; j++) 

if (omega[j] != -1) numer A = GF_table[ ((q-l-beta[i])*j + omega[j])%n ]; 

numer = dec_table[numer]; 

/* Correct the erroneous value */ 
if( numer != -1) 

r[beta[i]] A = GF_table[(n+numer+beta[il-den)%n]; 


} 

/* Change r back into the GF representation */ 
for(i=0; i<n; i++) r[i] = dec_table[r[i]]; 
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} /* end if (num_errors == deg_lambda) */ 

else 

{ 

decode_flag = 1 ; 

if(num_erasures == max_eras) counter[3] += 1; 

} 

} /* end if (num_erasures + num errors <= dmin) */ 

else 

{ 

decode_flag = 1 ; 

if(num_erasures == max_eras) counter[2] += 1 ; 


} 

} /* end if (degjambda <= 2*t) */ 

else 

{ 

decode_flag = 1 ; 

if(num_erasures == max_eras) counter[l] += 1; 

} 

} /* end if (error) */ 

else 

{ 

erasure_counter[num_erasures] += 1 ; 
decode_flag - 0; 
counter[0] += 1 ; 

} 


} /* end if decode flag */ 

} /* end for num_erasures = 0 to max_eras */ 


if (decode__flag) 

{ 

/* Incorrect decoding. */ 

counter[5] += 1 ; 

for(j=0; j<k; j++) 

{ 

for(l=0; l<m; 1++) 

data[ii*k*m + j*m + 1] = data[ii*n*m + (n-k+j)*m + 1]; 

} 

} 

else 

{ 

counter[4] += 1 ; 

/* Correct decoding */ 

for(j=0 ; j<k; j++) 

{ 

if (r[n-k+j] == -1) sum = 0; 
else sum = GF_tabIe[r[n-k+j]]; 
for (1=0; l<m; 1++) 
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{ 

cl - pow(2, m-1-1); 
if (sum >= cl) 

{ 

data [ii*k*m + j*m + 1] = 1; 
sum -= c 1 ; 

} 

else data [ii*k*m + j*m + 1] = 0; 

} 

} 

} 


free(GF_table); 

free(dec_table); 

free(reliability); 

free(min); 

free(beta); 

free(r); 

free(U); 

free(S); 

free(T); 

free(B); 

free(tmp); 

free(lambda); 

free(omega); 

} 


void rs_decode_erasure_method_2(int *data, int *v2, int ^counter, int n, int k, int q, int m, int num, 
int *erasure_counter, int I, int *decoded_word, 
int *frame_failure, int *error_stat 1 , int *error_stat2) 

/* Reed Solomon decoding is performed using erasure method 2. In this method, the deinterleaving */ 
/* frame is reconstructed. Decoding of each RSW is initially attempted using erasure method 1. */ 

/* The flags for each undecoded RSW are updated using information provided by the successfully */ 

/* decoded RSW in the frame. */ 

[ 

int GF_poly[l 1] = {0,0,0,6,12,20,48,72,184,272,576}; 
int i, ii, j, 1, 11, t, rr, sum, sum2, L, cl, degjambda, delta_r; 
int error, decode_flag, num_errors, numer, den, I_flag; 
int num_erasures, temp, max_eras, max_erasl, cc, failure, ne; 
int *r, *r2, *U, *S, *lambda, *omega, *B, *T, *tmp, *beta; 

struct mat **RSM; 

int jj, kk, reli, GEBE, offset, GE_count, r_count, max_iter, id, idl; 
int *GF_table, *dec_table; 

GF_table = ivector(n+l); 
dec_table = ivector(n+l); 

t = (n-k)/2; 
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beta = ivector(n); 
r = ivector(n); 
r2 = ivector(n); 

U = ivector(2*t); 

S = ivector(n); 

T = ivector(n); 

B = ivector(n); 
tmp = ivector(n); 
lambda = ivector(n); 
omega = ivector(2*t+l); 
counter = ivector(20); 

max_erasl = 16; 
max_iter = I; 

/* Create GF(2 A m) field */ 

make_GF_table (m, GF_table, dec_iable, GF_poly[m]); 

RSM = (struct mat **)calloc(n*m, sizeof(struct mat *)); 
for(i=0;i<n*m;i++) 

{ 

RSM[i] = (struct mat *)calloc(I, sizeof(struct mat)); 

} 


for (ii = 0; ii < num/I; ii ++) 

{ 


/* Enter the data into the matrix */ 

counter[13] ++; /* Total number of frames */ 

I_flag = 0; 

for(i=0; i<I; i++) 

{ 

id = 0; 

decoded_word[i] = 0; 
ne = 0; 

counter[0]++; /* number of RSW */ 

for(j = 0; j < n; j++) 

( 

RSM[j][i].r = 0; 

RSM[j][i].r2 = 0; 

RSM[J3[i3-reli = 500; 

RSM[j][i].flag = 0; 
sum = 0; 
sum2 = 0; 

for (1 = 0; 1 < m; 1++) 

{ 

if( RSM[j][i].reli > abs(data[ii*I*n*m + i*n*m + j*m +1])) 
RSM[j][i].reli = abs(data[ii*I*n*m + i^n^m + j*m +1]); 
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if(data[ii*I*n*m + i*n*m + j*m +1] <= 0) 
data[ii*I*n*m + i*n*m + j*m +1] = 0; 
else 

data[ii*I*n*m + i*n*m + j*m +1] = 1; 

, sum += data[ii*I*n*m + i*n*m + j*m +l]*pow(2, m-1-1); 
sum2 += v2[ii*I*n*m + i*n*m + j*m +l]*pow(2, m-1-1); 

} 

RSM[j][i].r = dec_table[sum]; 

RSM[j][i].r2 = dec_table[sum2]; 

++counter[ 1 1]; /* Total number of bytes */ 

if(sum != sum2) 

{ 

++ne; 

counter[12] ++; /* Number of byte errors */ 

} 


} 

if(ne >32) ne = 32; 

if(ne <=16) ne = 16; 
else I_flag = 1 ; 


++ error_statl[ne-16]; 

} 

if(I_flag) counter[14] ++; /* Number of frame errors */ 


/* If GEN_STAT - 1 , then the probability of flagging a BE given a GE */ 
/* and the probability of flagging a GE given a BE */ 

if(GEN_STAT) 

{ 

for(i=0; i<I; i++) 

{ 

for(j=0; j<n; j++) 

reli = RSM[j][i].reli; 
if(reli >= 255) continue; 

if(RSM[j][i].r == RSM[j][i].r2) GEBE = 2; 
else GEBE = 1; 

kk = j; 
offset = 0; 

/* loop forward looking for identical reliabilities */ 

for(l=0; Id; 1++) 

{ 
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if(i+l>=I && kk == n-1) break; 
if(i+I>=I) 


{ 

kk = j+1; 
offset = I; 

} 

if(RSM[kk][i+l-offset].reli = reli) 


if(GEBE == 1) 


counter[6]++; 

if(RSM[kk][i+l-offset].r == RSM[kk][i+l-offset].r2) 
counter[5]++; 
else counter[4]++; 


if(GEBE == 2) 


counter[9]++; 

if(RSM[kk][i+l-offset].r == RSM[kk][i+l-offset].r2) 
counter[7]++; 
else counter[8]++; 

} 

} 

else break; 

} 


/* loop backwards looking for identical identities */ 

kk = j; 
offset = 0; 

for(l=0; 1<I; 1++) 

{ 

if(i-l < 0 && kk == 0) break; 
if(i-l < 0) 

{ 

kk = j-1; 
offset = I; 

} 

if(RSM[kk][i-l+offset].reli == reli) 

{ 

if(GEBE = 1) 

{ 

counter[6]++; 

if(RSM[kk][i-l+offset].r == RSM[kk][i-l+offset].r2) 
counter[5]++; 
else counter[4]++; 


if(GEBE == 2) 

{ 

counter[9]++; 

if(RSM[kk][i-l+offset].r == RSM[kk][M+offset].r2) 
counter[7]++; 
else counter[8]++; 

} 
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} 

else break; 

} 

} /* end loop j = 0 to n */ 
} /* end loop i = 0 to I */ 

} /* end if(GEN_STAT) */ 


if(DEBUG_DUMP) 

{ 

printf("\nDecoded Word\n\n"); 

for(i=0; i<I; i++) printf("%d\t",decoded_word[i]); 

printf("\n"); 
for(i=0; i<n; i++) 

{ 

printf("\n%d ",i); 
for(j=0; jd; j++) 

{ 

if(RSM[i][j].r != RSM[i][j].r2) 
printf("(%d, GE)\t'\RSM[i][j].reli); 
else 

printf("(%d, BE)\t",RSM[ij[j].reli); 


} 

if(PAUSE) getchar(); 

} 


if(DEBUGl) 

{ 


printf("\nReceived word\n"); 

for(i=0; id; i++) 

{ 

printf("\n"); 
for(j=0; j<n; j++) 

{ 

printf(" %d",RSM[j][i].r); 

} 

} 


if(PAUSE) getchar(); 


printf("\nCorrect word\n"); 

for(i=0; id; i++) 

{ 

printf("\n"); 

for(j=0; j<n; j++) 

( 

printf(" %d",RSM[J][i].r2); 

} 
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} 


if(PAUSE) getcharQ; 


printf("\nReliability\n"); 

for(i=0; i<I; i++) 

{ 

printf("\n"); 
for(j=0; j<n; j++) 

{ 

printff %d'\RSM[j][i].reli); 

} 

} 

if(PAUSE) getchar(); 


printf("\nError flag\n"); 

for(i=0; i<I; i++) 

{ 

printf("\n"); 
for(j=0; j<n; j++) 

{ 

printff %d",RSM[j][i].flag); 

} 

} 

if(PAUSE) getchar(); 

} /* end if (DEBUG!) */ 

idl =0; 
id = 2; 

/* Attempt the iteritive decoding *1 

for(jj=0; jj<max_iter; jj++) 

{ 

/* Obtain the reliability info from the enibhoring codewords from */ 
/* the interleaver and/or SOVA */ 

if(jj==0) max_eras =16; 
else max_eras = max_eras 1 ; 

if(idl==id) break; 

id = idl; 

if(DEBUG) 

{ 

printf("\ndecoded word = "); 

for(i=0; id; i++) printf("%d ", decoded_word[i]); 
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if(id == I) break; 


for(i=0; i<I; i++) 

{ 


if(decoded_word[i] != 1) continue; 

for(j=0; j<n; j++) 

{ 

reii = RSM[j][i].re!i; 
if(reli >= 255) continue; 

GEBE = RSMIJ] [i] .flag; 

kk = j ; 
offset = 0; 

/* loop forward looking for identical reliabilities */ 

for(l=0; Id; I++) 

{ 

if(i+l>=I && kk == n-1) break; 
if(i+l>=I) 

{ 

kk = j+1; 
offset = I; 

} 

if(RSM[kk][i+l-offset].reli = reli) 

{ 

if(decoded_word[i+l-offset] == 0) 
RSM[kk][i+l-offset].flag = GEBE; 

} 

else break; 

} 

/* loop backwards looking for identical identities */ 

kk = j; 
offset = 0; 

for(l-0; Id; 1++) 

{ 

if(i-l < 0 && kk == 0) break; 
if(i-l < 0) 

{ 

kk = j-1; 
offset = I; 

} 

if(RSM[kk][i-l+offset].reli == reli) 

{ 

if(decoded_word[i-l+offset] == 0) 

RSM[kk] [i-l+offset] .flag = GEBE; 

} 


169 



else break; 


} 

} /* end loop j = 0 to n */ 

decoded_word[i] = 2; I* flag this word as already giving GEBE info */ 

} 

/* Done revising reliability info */ 

/* Start prioritising the reliability info, and then begin decoding */ 

for(ll = 0; 11<I; 11++) /* Find good erasures (GE) first */ 

{ 

GE_count = 0; 

if(DEBUG) 

{ 

printf("\ndecoded word = "); 

for(i=0; id; i++) printf("%d ", decoded_word[i]); 

} 

if(decoded_word[ll] != 0) continue; 

for(i=0; i<n; i++) 

{ 

if(RSM[i][ll] -flag == 1) 

{ 

U[GE_count] = i; 

GE_count ++; 

} 

} 

if(GE_count < max_eras) 

{ 

r_count = 0; 
for(i=0; i<n; i++) 

{ 

if(GE_count + r_count < max_eras) 

{ 

if(RSM[i][ll].flag == 0) 

{ 

U[GE_count + r_count] = i; 
r_count++; 

} 

} 

if(GE_count + r_count == max_eras) 

{ 

for(j=GE_count; j<max_eras; j++) 

{ 

for(l=GE_count; l<max_eras-l ; 1++) 

{ 

if(RSM[U[l]][ll].reli > RSM[U[l+l]][ll].reli) 

{ 
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temp = U[l] ; 
U[1] = U[1+1]; 
U[l+1] - temp; 


r_count ++; 


if(GE_count + r_count > max_eras) 

{ 

if(RSM[U[max_eras- l]][ll].reli > RSM[i][ll].reli 
&& RSM[i][ll].flag == 0) 

{ 

U[max_eras-1] = i; 

for(j=max_eras- 1 ; j>=GE_count+ 1 ; j— ) 

{ 

if(RSM[U[j]][ll].reli < RSM[U[j-l]][iI].reli) 

{ 

temp = Ujj] ; 

U[j] = U[j- 1 ]; 

U[j-1] = temp; 

} 

else break; 

} 

} 

} /* End if(GE_count + r_count > max_eras) */ 

} /* End loop i = 0 .. n */ 

} /* End loop if(GE_count < max_eras) */ 

/* Now we have reliabilities try to decode */ 


if(DEBUG) 

{ 

printf("\nMinimum reliabilities and positions'^"); 

for(i=0; i<max_eras; i++) printf("(%d, %d)",U[i], RSM[U[i]][ll].reii); 

printf("\n"); 

printf("\nnumber of GE = %d", GE_count); 
cc - 0; 

for(i=0; i<n; i++) 

{ 

if(RSM[i][ll].r != RSM[i][ll].r2) 

{ 

cc += 1 ; 

printf("\n%d error at %d Reliability = %d", cc, i, RSM[i][lI].reli); 
for(j=0; j<max_eras; j++) 

{ 

if(i == U[j]) printf(" min %d", j); 


} 
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I 

if(PAUSE) getchar(); 

} /* end if (DEBUG) */ 

for(i=0; i<n; i++) r[i] = RSM[i][ll].r; 

/* Begin decoding. Try with 0 erasures, and keep on adding 2 */ 

/* until decodes properly, or until maximum is reached */ 

/* The symbols with the minimum reliabilities are in positions */ 

/* U[0] ... U[num_erasures-1], These symbols will be erased. */ 

/* Set these symbols equal to 0 before computing the syndrome */ 

decode_flag = 1 ; 

for(num_erasures=0; num_erasures <= max_eras; num_erasures +=2) 

{ 


if(decode_flag) 

{ 

counter[15] ++; /* Number of decoding trials */ 
for(i=0; i<num_erasures; i++) r[U[i]] = -l; /* =0 */ 
error = 0; 

/* Compute the syndrome */ 
for(i=l; i<= n-k; i++) 

{ 

S[i] = 0; 

for(j=0; j<n; j++) 

if(rU3 !- -1) 

S[i] A = GF_table[ (r[j] + i*j) % n]; 
if (S[i] != 0) error = 1; /* If nonzero syndrome, there is an error*/ 

} 

S[0] = 0; 

/* Convert to GF representation */ 
for(i=l; i<=n-k; i++) S[i] = dec_table[S[i]]; 

if (error) /* If the syndrome is equal to zero, no decoding nessasary */ 

{ 

for(i=0; i < n; i++) 

{ 

lambda[i] = 0; 

B[i] = -1; 

T[i] = 0; 

} 

lambda[0] = GF_table[0], /* =1 */ 

L = 0; 
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degjambda = 0; 

B [0] = 0; /*=!*/ 


for(rr= 1 ; rr <= 2*t; rr++) 

{ 

if (rr <= num_erasures 1 ) 

{ 

for (i=l; i <= degjambda+l; i++) 

{ 

if (lambda[i-l] != 0) 

tmp[i] = GF_table[(U[rr-l]+dec_table[lambda[i-l]])%n]; 
else tmp[i] = 0; 

} 

for 0=1; j <= degjambda+l; j++) lambda[j] A = tmp[j]; 
degjambda++; 

for(j=0; j<= degjambda; j++) B[j] = decjable[lambda[j]]; 
L ++; 

} 

else 

{ 

/* Compute the discrepancy */ 
delta_r = 0; 
for(j=0; j<=rr; j++) 
if (lambdajj] != 0 && S[rr-j] !=-l) 
delta_r A = GF_table[(dec_table[lambda(j]]+S[rr-j])%n]; 
delta_r = decjable[delta_r]; 

if (delta_r != -1) 

{ 

/* T(x) <— Iambda(x) + delta_r * x * B(x) */ 

for (j= 1 ; j <- deg_lambda + 1 ; j++) 
if (B(j-1] != -1) 

T[j] = lambdafj] A GF_table[ (delta_r + B[j-1]) % n]; 
T[0] = lambda[0]; 

++deg_lambda; 

if (2*L <= rr + num_erasures-l) 

{ 

L = rr - L + num_erasures; 

/* B(x) < — B(x) / delta_r */ 
for(i=0; j<= degjambda; j++) 
if( lambdaQ] != 0) 

B[j] = (n-delta_r+decjable[lambda[j]])%n; 
else B[j] = -1; 

/* lambda(x) < — T(x) */ 

for(j=0; j<=deg Jambda+ 1 ; j++) lambda[j] = T[jJ; 


} 

else 

{ 

/* lambda(x) < — T(x) */ 
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for(j=0; j<n-k; j++) lambda!]] = T[j]; 


/* B(x) < x * B(x) */ 

tmp[0] = - 1 ; 

for(j= 1 ; j<=n-k; j++) tmp[j] = B[j-1]; 
for(j=0; j<=n-k; j++) B[j] = tmp[j]; 

} 


else 

{ 

/* B(x) < x * B(x) */ 

tmp[0] = - 1 ; 

for(j= 1 ; j<=n-k; j++) tmp[j] = B[j-1]; 
■ for(J=0; j<=n-k; j++) B[j] = tmp[j]; 

} 

} 

} 


/* Change lambda(x) to GF representation */ 
for(i=0; i<n; i++) lambda[i] = dec_table[lambda[i]]; 

/* Compute the degree of lambda(x) */ 
deg_lambda = n; 
for(i=n- 1 ; i>=0; i--) 

if (lambda[i] != -1 && deg_lambda == n) deg_lambda = i; 

if (degjambda <= 2*t) /* Below the capacity of the code */ 

{ 


/* Comupte omega(X) = [1+S(X)] * lambda(X) */ 


for(i=0; i<=n-k; i++) omega[i] = 0; 
for (i=0; i<=n-k; i++) 

{ 

for 0=0; j<=n-k; j++) 


if ((i+j) >= n-k+1) continue; 

if (S[i] !- -1 && lambda!]] != -1) 

omega[i+j] A = GF_table[ (S[i] + lambda!]]) % n ]; 

} 


/* Convert omega(x) to GF representation */ 

for (i=0; i<=n-k; i++) omega[i] = dec_table[omega[i]]; 

/* Find the roots of lambda(X). The inverses of the roots gives us the */ 
/* location of the errors. */ 

num_errors = 0; 
for (i=0; i<q- 1 ; i++) 

{ 


sum = 0; 

for (j=0; j<2*t; j++) 
if (lambda[j] !=-l) 
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sum A = GF_table[ (lambda[j] + i*j) % (q~l) ]; 
if (sum == 0) 

{ 

beta[num_errors] = (n-i)%n; 
num_errors++; 


} 

if( (2*(num_errors-num_erasures) + num_erasures) <= 2*t) 

{ 

if (num_errors == deg_lambda) 

{ 


/* Correct Decoding, record the umber of erasures required */ 
erasure_counter[num_erasures] += 1 ; 
decode_flag = 0; 

/* Convert r to base 1 0 representation */ 

for(i=0; i<n; i++) 
if (r[i] != -1) r[i] = GF_table[r[iJ]; 
else r[i] = 0; 

/* Calculate the error values and correct the received vector */ 

for (i=0; i<num„errors; i++) 

{ 

/* Calculate the denominator */ 
den = 0; 

for (j=0; j<=2*t; j++) 


if (j%2 == 0) continue; 
if (lambda[j] != -1) 

{ 

cl = GF_table[ ( (q-l-beta[i])*(j-l) + lambdaQ]) % n ]; 
den A = c 1 ; 

} 


den = dec_table[den]; 

/* Calculate the numerator */ 

numer = 0; 

for (j=0; j<=2*t; j++) 

if (omega[j] != -1) numer A = GF_table[ ((q-l-beta[i])*j + omega[j])%n ]; 

numer = dec_table[ numer]; 

/* Correct the erroneous value */ 
if(numer != -1) 

r[beta[i]] A = GF_table[(n+numer+beta[i]-den)%n]; 
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} 


/* Change r back into the GF representation */ 

for(i=0; i<n; i++) r[i] = dec_table[r[i]]; 

} /* end if (num_errors == degjambda) */ 

else 

{ 

decode_flag = 1 ; 

} 

} /* end if (num_erasures + num errors <= dmin) */ 

else 

{ 

decode_flag = 1 ; 

} 

} /* end if (deg_lambda <= 2*t) */ 

else 

{ 

decode_flag = 1 ; 

} 

} /* end if (error) */ 

else 

{ 

erasure_counter[num_erasures] += 1 ; 
decode_flag = 0; 

} 

} /* end if decode flag */ 


} /* end for num_erasures = 0 to max_eras */ 


if (!decode_flag) 

{ 

/* Correct decoding. */ 


decoded_word[ll] = 1; 
idl +=1; 

for(i=0; i<n; i++) 

{ 

if(RSM[i][ll].r != RSM[i][ll].r2) RSM[i][ll].flag = 1; /* GE */ 
else RSMfi] [11] .flag = 2; 

} 

if(DEBUG) 

{ 

printf("\nWord %d decoded correctly", 11); 
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if(PAUSE) getchar(); 


} 

failure = 0; 
for(i=0; i<n; i++) 

{ 

RSM[i][ll].r = r[i]; 

if(RSM[i][ll].r != RSM[i][ll].r2) failure = 1; 

} 


if(failure) counter[10]++; 


} 

} /* End 11 = 0 ... I */ 

} /* End jj = 0 .. max_iter */ 


ne = 0; 

for(jj=0; jj<I; jj++) if(decoded_word[jj]!=0) ++ne; 
++frame_failure[ne] ; 

/* Send the result to the output */ 


for(jj=0; jj<I; jj++) 


if(decoded_word[jj] == 0) 

{ 

++counter[2]; /* Number of incorrect RSW */ 
ne = 0; 

for(j=0; j<n; j++) 

if(RSM[j][jj].r != RSM[j][jj].r2) ne++; 
++error_stat2[ne- i 6] ; 

} 

else ++counter[l]; /* Number of correct RSW */ 

for(j=0 ; j<k; j++) 

{ 

if (RSM[n-k+j][jj].r == -1) sum = 0; 
else sum = GF_table[RSM[n-k+j][jj].r]; 
for (1=0; l<m; 1++) 

{ 

cl = pow(2, m-1-1); 
if (sum >= cl) 

{ 

data [ii*I*k*m + +jj*m*k + j*m + 1] = 1; 
sum -= c 1 ; 

} 

else data [ii*I*k*m + +jj*m*k + j*m + 1] = 0; 


} 
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for(i=0;i<n*m;i++) 

free(*RSM); 

free(RSM); 

free(beta); 

free(r); 

free(r2); 

free(U); 

free(S); 

free(T); 

free(B); 

free(tmp); 

free(lambda); 

free(omega); 

free(counter); 

} 


void rs_decode_erasure_method_3(int *data, int *v2, int *counter, int n, int k, int q, 
int m, int num, int *erasure_counter, int I, int *decoded_word, 
int *frame_failure, int *error_statl, int *error_stat2) 

{ 

int GF_poly[ll] = {0,0,0,6,12,20,48,72,184,272,576}; 

int i, ii, j, 1, 11, t, it, sum, sum2, L, cl, degjambda, delta_r; 

int error, decode_flag, num_errors, numer, den, I__flag; 

int num_erasures, temp, max_eras, max_erasl, cc, failure, ne; 

int *r, *r2, *U, *S, *lambda, *omega, *B, *T, *tmp, *beta; 

int jj, kk, reli, GEBE, offset, GE_count, r_count, max_iter, id, idl; 

int redecode_flag, new_info_flag; 

int *GF_table, *dec_table; 

struct mat **RSM; 

GF_table = ivector(n+l); 
dec_table = ivector(n+l); 

t = (n-k)/2; 

beta = ivector(n); 
r = ivector(n); 
r2 = ivector(n); 

U = ivector(2*t); 

S = ivector(n); 

T = ivector(n); 

B = ivector(n); 
tmp = ivector(n); 
lambda = ivector(n); 
omega = ivector(2*t+l); 


max_erasl = 16; 
max_iter = I; 
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/* Create GF(2 A m) field */ 


make_GF_table (m, GF_table, dec_table, GF_poly[m]); 


for (ii = 0; ii < num/I; ii ++) 

{ 


/* Enter the data into the matrix */ 

counter[13] ++; /* Total number of frames */ 

I_flag = 0; 

for(i=0; i<I; i++) 

{ 

id = 0; 

decoded_word[i] = 0; 
ne = 0; 

counter[0]++; /* number of RSW */ 

for(j = 0; j < n; j++) 

{ 

RSM[j][i].r = 0; 

RSM[j][i].r2 = 0; 

RSM[j][i].reli = 500; 

RSM[j][i].flag = 0; 
sum = 0; 
sum2 = 0; 

for (1 = 0; 1 < m; 1++) 

{ 

if( RSM[j][i].reli > abs(data[ii*I*n*m + i*n*m + j*m +1])) 
RSM[j][i].reli = abs(data[ii*I*n*m + i*n*m + j*m +1]); 

if(data[ii*I*n*m + i*n*m + j*m +1] <= 0) 
data[ii*I*n*m + i*n*m + j*m +1] = 0; 
else 

data[ii*I*n*m + i*n ;f: m + j*m +1] = 1; 

sum += data[ii*I*n*m + i*n*m + j*m +l]*pow(2, m-1-1); 
sum2 += v2[ii*I s|t n*m + i*n*m + j*m +l]*pow(2, m-1-1); 

} 

RSM[j][i].r = dec_table[sum]; 

RSM[j][i].r2 = dec_table[sum2]; 

++counter[llj; /* Total number of bytes *1 

if(sum != sum2) 

{ 

++ne; 

counter[12] ++; /* Number of byte errors */ 

} 


} 

if(ne >trunc_length) ne = 32; 
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if(ne <=16) ne = 16; 
else I_flag = 1 ; 


++ error_statl[ne-16]; 

} 


if(I_flag) counter[14] ++; /* Number of frame errors */ 


/* If GEN_STAT = 1, then the probability of flagging a BE given a GE */ 
/* and the probability of flagging a GE given a BE */ 

if(GEN_STAT) 

{ 

for(i=0; i<I; i++) 

{ 

for(j=0; j<n; j++) 

{ 

reli = RSM(j][i].reli; 
if(reli >= 255) continue; 

if(RSM[j][i].r == RSM[j][i].r2) GEBE = 2; 
else GEBE = 1; 

kk = j; 
offset = 0; 

/* loop forward looking for identical reliabilities */ 

for(!=0; Id; 1++) 

{ 

if(i+l>=I && kk == n-1) break; 
if(i+l>=I) 

{ 

kk = j+1; 
offset = I; 

} 

if(RSM[kk][i+l-offset].reli == reli) 


if(GEBE == 1) 

{ 

counter[61++; 

if(RSM[kk] [i+l-offset] .r == RSM[kk][i+l-offset].r2) 
counter[5]++; 
else counter[4]++; 

} 

if(GEBE == 2) 


counter[9]++; 

if(RSM[kk][i+l-offset].r == RSMfkk] [i+l-offset]. r2) 
counter[7]++; 
else counter[8]++; 

} 
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} 

else break; 

} 

/* loop backwards looking for identical identities */ 

kk = j; 
offset = 0; 

for(l=0; 1<I; 1++) 

{ 

if(i-l < 0 && kk == 0) break; 
if(i-l < 0) 

{ 

kk = j-1; 
offset = I; 

} 

if(RSM[kk][i-l+offset].reli == reli) 

{ 

if(GEBE — 1) 

{ 

counter[6]++; 

if(RSM[kk][i-l+offset].r == RSM[kk][i-l+offset].r2) 
counter[5]++; 
else counter[4]++; 


if(GEBE == 2) 

{ 

counter[9]++; 

if(RSM[kk][i-l+offset].r == RSM[kk][i-l+offset].r2) 
counter[7]++; 
else counter[8]++; 


else break; 

} 

} /* end loop j = 0 to n */ 
} /* end loop i = 0 to I */ 

} /* end if(GEN_STAT) */ 


if(DEBUG_DUMP) 

{ 

printf("\nDecoded Word\n\n"); 

for(i=0; id; i++) printf(' , %d\t",decoded_wordn]); 

printfCNn"); 
for(i=0; i<n; i++) 

{ 

printf("\n%d ”,i); 
for(j=0; jd; j++) 

{ 

if(RSM[i][j].r != RSM[i](j].r2) 
printf("(%d, GE)\t",RSM[i][j].reli); 
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else 

printf("(%d, BE)\t",RSM[i][j].reli); 

} 

} 

if(PAUSE) getchar(); 

} 


if(DEBUGl) 

{ 

printf("\nReceived word\n"); 

for(i=0; i<I; i++) 

{ 

printf("\n"); 
for(j=0; j<n; j++) 

{ 

printfC’ %d",RSM[j][i].r); 

} 

} 

if(PAUSE) getchar(); 


printf("\nCorrect word\n"); 

for(i=0; i<I; i++) 

{ 

printf("\n"); 
for(j=0; j<n; j++) 

{ 

printf(" %d",RSM(j][i].r2); 

} 

} 

if(PAUSE) getchar(); 


printf("\nReliability\n"); 

for(i=0; i<I; i++) 

{ 

printf("\n"); 
for(j=0; j<n; j++) 

{ 

printfC %d",RSM[j][i].reli); 

} 

} 

if(PAUSE) getchar(); 


printf("\nError flag\n”); 
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for(i=0; i<I; i++) 


( 

printf("\n"); 
for(j=0; j<n; j++) 

{ 

printff %d",RSM[j][i].flag); 


} 

if(PAUSE) getchar(); 

} /* end if (DEBUG 1)*/ 

idl =0; 
id = 2; 

/* Attempt the iteritive decoding */ 

f°r(jj=0; jj<max_iter; jj++) 

{ 


/* Obtain the reliability info from the enibhoring codewords from */ 
/* the interleaver and/or SOYA */ 


if(jj--0) max_eras =16; 
else max_eras = max_eras l ; 

if(idl==id) break; 

id = id 1 ; 

if(DEBUG) 

{ 

printf("\ndecoded word = "); 

for(i=0; id; i++) printf("%d ", decoded_word[i]); 

} 


if(id == I) break; 

if(jj==0) redecode_flag = 1 ; 
else redecode_flag = 0; 

for(i=I-l; i>=l; i— ) 

if(decoded_word[i] !=0 && decoded_word[i-l] == 0) 
redecode_flag = 1 ; 

if(!redecode_flag) break; 

/* Start prioritising the reliability info, and then begin decoding */ 

for(ll = 0; lid; 11++) /* Find good erasures (GE) first */ 

{ 

GE_count = 0; 

if(DEBUG) 

{ 
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printf("\ndecoded word = "); 

for(i=0; i<I; i++) printf("%d ", decoded_word[i]); 

} 


if(decoded_word[ll] != 0) continue; 

for(i=0; i<n; i++) 

{ 

if(RSM[i][ll] .flag — 1) 

{ 

U[GE_count] = i; 

GE_count ++; 

} 

} 

if(GE_count < max_eras) 

{ 

r_count = 0; 
for(i=0; i<n; i++) 

{ 

if(GE_count + r_count < max_eras) 

{ 

if(RSM[i][ll] .flag = 0) 

{ 

U[GE_count + r_count] = i; 
r_count++; 

} 

} 


if(GE_count + r_count == max_eras) 

{ 

for(j=GE_count; j<max_eras; j++) 

{ 

for(l=GE_count; l<max_eras- 1 ; 1++) 

{ 

if(RSM[U[l]][ll].reli > RSM[U[l+l]][U].reli) 

{ 

temp = U[l]; 

U[1] = U[1+1]; 

U[l+1] = temp; 


} 

} 

r_count ++; 

} 

if(GE_count + r_count > max_eras) 

{ 

if(RSM[U[max_eras-l]][ll].reli > RSM[i][ll].reli 
&& RSM[i][il].flag == 0) 

{ 


U[max_eras-1] = i; 
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for(j=max_eras- 1 ; j>=GE_count+ 1 ; j— ) 

{ 

if(RSM[U(j]][ll].reli < RSM[U[j-l]][ll].reli) 

{ 

temp = U[j]; 

U(j] = UO-l]; 

U[j-13 = temp; 

} 

else break; 

} 

} /* End if(GE_count + r_count > max_eras) */ 

} /* End loop i = 0 .. n *1 
} /* End loop if(GE_count < max_eras) */ 

/* Now we have reliabilities try to decode */ 


if(DEBUG) 

{ 

printf("Number of trials = %d\n\nMinimum reliabilities and positions\n", counter[15]); 
for(i=0; i<max_eras; i++) printf("(%d, %d)",U[i], RSM[U[i]][ll].reIi); 
printf("\n"); 

printf("\nnumber of GE = %d'\ GE_count); 
cc = 0; 

for(i=0; i<n; i++) 

{ 

if(RSM[i][ll].r != RSM[i][ll].r2) 

{ 

cc += 1 ; 

printf("\n%d error at %d Reliability = %d", cc, i, RSM[i][ll].reli); 
for(j=0; j<max_eras; j++) 

{ 

if(i == U[j]) printf(" min %d", j); 

} 

if (RSM[i][ll].flag == 1) printf(" GE"); 


if(PAUSE) getchar(); 

} /* end if (DEBUG) */ 

for(i=0; i<n; i++) r[i] = RSM[i][lI].r; 

decode_flag = 1 ; 

if(GE_count >= max_erasl) max_eras = GE_count; 
else max_eras = max_eras 1 ; 

for(num_erasures=GE_count; num_erasures <= max_eras; num_erasures +=2) 

{ 
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if(decode_flag) 

{ 

counter[15] ++; /* Number of decoding trials */ 
for(i=0; i<num_erasures; i++) r[U[i]] = -l; /* =0 */ 
error = 0; 

/* Compute the syndrome */ 
for(i=l; i<= n-k; i++) 

{ 

S[i] = 0; 

for(j=0; j<n; j++) 
if(r[j]!=-l) 

S[i] A = GF_table[ (r[j] + i*j) % n]; 
if (S[i] != 0) error = 1 ; /* If nonzero syndrome, there is an error*/ 

} 

S[0] - 0; 

/* Convert to GF representation */ 
for(i= 1 ; i<=n-k; i++) S[i] = dec_table[S[i]]; 

if (error) /* If the syndrome is equal to zero, no decoding nessasary */ 

{ 

for(i=0; i < n; i++) 

{ 

lambdafi] = 0; 

B[i] = -1; 

T[i] = 0; 

} 


lambdafO] = GF„table[0]; /*={*/ 

L = 0; 

deg_lambda = 0; 

B[0] = 0; /* = 1 '*/ 

for(rr= I ; rr <= 2*t; rr++) 

{ 

if (rr <= num_erasures) 

{ 

for (i=l; i <= deg_lambda+l ; i++) 

{ 

if (lambda[i-lj !=0) 

tmp[i] = GF_table[(U[rr-l]+dec_table[lambda[i-l]])%n]; 
else tmp[i] = 0; 

} 

for (j=l ; j <= deg_lambda+ 1 ; j++) lambda(j] A = tmp[j]; 
deg_lambda++; 

for(j=0; j<= deg_lambda; j++) B[j] = dec_table[lambda[j]]; 
L ++; 

} 


else 

{ 
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/* Compute the discrepancy */ 
delta_r = 0; 
for(j=0; j<=rr; j++) 
if(lambda[j] != 0 && S[rr-j] !=-l) 
delta_r A = GF_table[(dec_table[lambda[j]]+S[rr-j])%n]; 
delta_r = dec_table[delta_r]; 

if (delta_r != -1) 

{ 

/* T(x) < — lambda(x) + delta_r * x * B(x) */ 

for (j= 1 ; j <= degjambda + 1 ; j++) 
if (B[j- 1 ] != -1) 

T[j] = lambda[j] A GF_table[ (delta_r + B [j- 1 ]) % n] ; 
T[0] = lambda[0]; 

++deg_lambda; 

if (2*L <= rr + num_erasures-l) 

{ 

L = rr - L + num_erasures; 

/* B(x) <— B(x) / delta_r */ 
for(j=0; j<= deg_lambda; j++) 
if( lambda[j] != 0) 

BQ] = (n-delta_r+dec_table[lambda[]]])%n; 
else B [j] = -1; 

/* lambda(x) < — T(x) */ 

for(j=0; j<=deg_lambda+ 1 ; j++) lambdafj] = TJJ]; 


} 

else 

{ 

/* lambda(x) < — T(x) */ 
for(j=0; j<n-k; j++) lambda(j] = T[j]; 

/* B(x) < x * B(x) */ 

tmp[0] = -1; 

for(j=l; j<=n-k; j++) tmp(j] = B[j-1]; 
for(j=0; j<=n-k; j++) B[j] = tmp[j]; 

} 


else 

{ 

/* B(x) < x * B(x) */ 

tmp[0] = -1; 

for(j= 1 ; j<=n-k; j ++) tmp[j] = B[j-1]; 
for(j=0; j<=n-k; j++) B[j] = tmpOl; 

} 

} 

) 


/* Change lambda(x) to GF representation */ 
for(i=0; i<n; i++) lambda[i] = dec__table[lambda[i]]; 
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/* Compute the degree of lambda(x) */ 
deg_lambda = n; 
for(i=n-l; i>=0; i--) 

if (lambda[i] != -1 && degjambda == n) deg_lambda = i; 

if (deg_lambda <= 2*t) /* Below the capacity of the code */ 

{ 


/* Comupte omega(X) = [l+S(X)j * lambda(X) */ 

for(i=0; i<=n-k; i++) omega[i] = 0; 
for (i=0; i<=n-k; i++) 

{ 

for (j— 0; j<=n-k; j++) 


if ((i-fj) >= n-k+1) continue; 

if (S[i] != -1 && lambdaQ] !=-l) 

omega[i+j] A = GF_table[ (S[i] + lambda[j]) % n ]; 

} 


/* Convert omega(x) to GF representation */ 

for (i=0; i<=n-k; i++) omegafi] = dec_table[omega[i]]; 

/* Find the roots of lambda(X). The inverses of the roots gives us the */ 
/* location of the errors. */ 

num_errors = 0; 
for (i=0; i<q- 1 ; i++) 

{ 

sum = 0; 

for (j=0; j<2*t; j++) 
if (lambda[j] != -1) 

sum A = GF_table[ (lambda[j] + i*j) % (q-1) ]; 
if (sum == 0) 


beta[num_errors] = (n-i)%n; 
num_errors++; 


if( (2*(num_errors-num_erasures) + num_erasures) <= 2*t) 

{ 

if (num_errors == degjambda) 

{ 

/* Correct Decoding, record the umber of erasures required */ 
erasure_counter[num_erasures] += 1; 
decode_flag = 0; 

/* Convert r to base 10 representation */ 
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for(i=0; i<n; i++) 
if (r[i] !=-l) r[i] = GF_table[r[i]]; 
else r[i] = 0; 


/* Calculate the error values and correct the received vector */ 

for (i=0; i<num_errors; i++) 

{ 

/* Calculate the denominator */ 
den = 0; 

for (j=0; j<=2*t; j++) 

{ 

if (j%2 == 0) continue; 
if (lambda[j] != -1) 

{ 

cl = GF_table[ ( (q-l-beta[i])*(j-l) + lambda[j]) % n ]; 
den A = c 1 ; 

} ■ 

} 

den = dec_table[den]; 

/* Calculate the numerator */ 

numer = 0; 

for (j=0; j<=2*t; j++) 

if (omega[j] != -1) numer A = GF_table[ ((q-l-beta[i])*j + omega[j])%n ]; 

numer = dec_table[numer]; 

/* Correct the erroneous value */ 
if(numer != -1) 

r[beta[i]] A = GF_table[(n+numer+beta[i]-den)%n]; 

} 

/* Change r back into the GF representation */ 

for(i=0; i<n; i++) r[i] = dec_table[r[i]]; 

} /* end if (num_errors == deg_lambda) */ 

else 

{ 

decode_flag = 1 ; 

} 

} /* end if (num_erasures + num errors <= dmin) */ 

else 

{ 

decode_flag = 1 ; 

} 

} /* end if (degjambda <= 2*t) */ 

else 

{ 
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decode_flag = 1 ; 

} 

} /* end if (error) */ 

else 

{ 

erasure_counter[num_erasures] += l; 
decode_fiag = 0; 

} 

} /* end if decode flag */ 


} /* end for num_erasures = 0 to max_eras */ 


if (!decode_flag) 

{ 

/* Correct decoding. */ 


decoded_word[ll] = 1; 
idl +=1; 


for(i=0; i<n; i++) 

{ 

if(RSM[i][llj.r != RSM[i][ll].r2) RSM[i][llj.flag = 1; /* GE */ 
else RSM[i][ll].flag = 2; 

} 


if(DEBUG) 

{ 

printf("\nWord %d decoded correctly", 11); 
if(PAUSE) getchar(); 

} 


failure = 0; 
for(i=0; i<n; i++) 

{ 

RSM[i][ll].r = r[i]; 

if(RSM[i][ll].r != RSM[i][ll].r2) failure = 1; 

} 


if(failure) counter[10]++; 
/* Begin updating the flags */ 


for(j=0; j<n; j++) 

{ 

reli = RSM[j][Il].reli; 
if(reli >= 255) continue; 
GEBE = RSM[j][ll].flag; 
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kk = j ; 
offset = 0; 

/* loop forward looking for identical reliabilities */ 

for(l=0; I<I; 1++) 

{ 

if(ll+l>=I && kk == n-1) break; 
if(ll+l>=I) 

{ 

kk = j+1; 
offset = I; 

} 

if(RSM[kk] [11+1-offset], reli == reli) 

{ 

if(decoded_word[ll+l-offset] == 0) 

RSM[kk] [11+1-offset], flag = GEBE; 

} 

else break; 

} 

/* loop backwards looking for identical identities */ 

kk = j; 
offset = 0; 

for(l=0; 1<I; 1++) 

{ 

if(ll-l < 0 && kk == 0) break; 
if(ll-l < 0) 

{ 

kk = j- 1 ; 
offset = I; 

} 

if(RSM[kk][U-l+offset].reli == reli) 

{ 

if(decoded_word[ll-l+offset] == 0) 

RSM[kk][ll-l+offset].flag = GEBE; 

} 

else break; 

} 

} /* end loop j = 0 to n */ 

decoded_word[i] = 2; /* flag this word as already giving GEBE info */ 


/* Done revising reliability info */ 


} /* End 11 = 0 ... I */ 

} /* End jj = 0 .. max_iter */ 
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ne = 0; 

for(jj=0; jj<I; jj++) if(decoded_word[jj]!=0) ++ne; 
++frame_failure[ne] ; 

/* Send the result to the output */ 


for(jj=0; jj<I; jj++) 

{ 

if(decoded_word[jj] == 0) 

{ 

++counter[2]; /* Number of incorrect RSW */ 
ne = 0; 

for(j=0; j<n; j++) 

if(RSM[jj|jjj.r != RSM[j][jj].r2) ne++; 
++error_stat2[ne- 1 6] ; 

} 

else ++counter[l]; /* Number of correct RSW *1 

for(j=0 ; j<k; j++) 

{ 

if (RSM[n-k+j][jj].r == -1) sum = 0; 
else sum = GF„table[RSM[n-k+j][jj3.r]; 
for (1=0; l<m; 1++) 

{ 

cl = pow(2, m-1-1); 
if (sum >= cl) 

{ 

data [ii*I*k*m + +jj*m*k + j*m + 1] = 1; 
sum -= cl; 

} 

else data [ii*I*k*m + +jj*m*k + j*m + 1] = 0; 



} 

for(i=0;i<n*m;i++) 

ffee(*RSM); 

free(RSM); 

free(beta); 

free(r); 

free(r2); 

free(U); 

free(S); 

free(T); 

free(B); 

free(tmp); 

free(lambda); 

free( omega); 

} 
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B.8 Main program for real system using erasure Method 1 


# include <stdio.h> 

# include <math.h> 

# include <stdlib.h> 

# include <stddef.h> 

# include "mmt31.h" 

# define L_DEBUG 0 

# define PAUSE 0 

# define Nss 16 /* Number of samples per symbol */ 

# define Rb 1 /* The bit rate */ 

# define Tb(1.0/Rb) 

# define Eb 1 .0 /* The Energy per bit */ 

# define Pi 3.14159265359 

/* Defining the random number generator constants */ 


# define IM1 2147483563 

# define IM2 2147483399 

# define AM ( 1 .0/IM 1 ) 

# define IMM1 (INI 1- 1 ) 

# define IA1 40014 

# define IA2 40692 

# define IQ1 53668 

# define IQ2 52774 

# define IR1 12211 

# define IR2 3791 

# define NTAB 32 

# define NDIV (1+IMM1/NTAB) 

# define EPS 1.2e-7 

# define RNMX(l.O-EPS) 


void main(int argc, char *argv[]) 

{ 

int *v, *u, i, ii, errors, Ml, M2, Q; 
long *idum; 
long idum_cell; 

long num_samples, num_bits, N, N_total; 
int rs_m, rs_t, rs_n, rs_k, rs_q, rs_num; 
int conv_k, conv_n, conv_m, conv_num; 
double *data, *Tx, *Rx, Pe; 
double Psig, EbNo, gain; 
int max_repeat, I, rows, I_num, a; 

double *xl, *x2; 
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FILE *inpfO, *inpfl, *inpf2, *inpf3, *outfO; 

int **soft_metric, *erasure__counter, *rs_counter; 
int num_erasures; 
char cFilename[80]; 

idum=&idum_ceil; 

EbNo=atof(arg v [ 1 ] ) ; 

N = atol(argv[2]); 

I = atol(argv[3]); 
max_repeat = atol(argv[4]); 


sprintf(cFilename,"New_Real_Ml_erasure„%2.2fEbNo_%dI_32delta_8soft.output",EbNo,I); 

inpfO = fopen("ml.dat","r"); 
inpfl = fopen("m2.dat","r"); 
inpf2 = fopen("tx.dat","r"); 
inpf3 = fopen(''rx.dat","r"); 


/* obtain the number of single sided filter */ 

/* coefficents for Tx (Ml) and Rx (M2) */ 

fscanf(inpfO, "%d", &M1); 
fscanf(inpfl, "%d", &M2); 

Tx = dvector(2*Ml); 

Rx = dvector(2*M2); 

xl = dvector(2*Ml+l); 
x2 = dvector(2*M2+l); 

/* Read in the filter coefficients from data files */ 

for(i=0; i<2*Ml; i++) 
fscanf(inpf2, &Tx[i]); 

for(i=0; i<2*M2; i++) 
fscanf(inpf3, "%lf', &Rx[i]); 


/* The Reed Solomon coding parameters */ 

rs_m = 8; 
rs_t = 16; 

rs_q = pow(2, rs_m); 
rs_n = rs_q- 1 ; 
rs_k = rs_n - 2*rs_t; 

/* Interleaving parameters */ 


rows = rs_m*rs_n; 
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/* The Convolutional coding parameters */ 

conv_k= 1; 
conv_n = 2; 
conv_m - 6; 

Q = B; /* Number of soft decision levels */ 

rs_num = N/(rsJk*rsjm); 

if( N % (rs_k*rs_m) != 0) 

{ 

++rs_num; 

N = rs_num*rs_k*rs_m; 

} 

I_num = rs_num/I; 


if(rs_num%I != 0) 

{ 

++ I_num; 
rs_num = I_num*I; 

N = rs_num*rs_k*rs_m; 

} 


gain = (double)(conv_k*rs_k)/(conv_n*rs_n); 
a = 10.0; 

conv_num = (rs_num*rs_n*rs_m + 50)/conv_k; 
num_bits = conv_n*(conv_num + 100); 
num_samples - num_bits*Nss; 
errors = 0; 

erasure_counter = ivector(2*rs_t); 
rs_counter = ivector(rs_t); 

soft_metric = int_matrix_2d(2, Q); 

for(i=0; i<6; i++) rs_counter[i] = 0; 

N_totaI - 0; 

data = dvector(num_samples+3*M 1 ); 


v = ivector(num_bks); 
u = ivector(N); 


*idum = -10; 
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for(ii=0; ii<max_repeat; ii++) 

{ 

for (i=0; i<num_bits; i++) v[i] = 0; 


for (i=0; i<N; i++) 

{ 

v[i] = bitgen(idum); 
u[i] = v[i] ; 

} 

rs_encode(v, rs_n, rs_k, rs_m, rs_num); 
interleave(v, rows, rs_n, rs_m, I, I„num); 
conv_encode(v, conv_n, conv_k, conv„m, conv_num); 
modulate(data, v, num_bits); 


filter(data, Ml, Tx, xl, num_samples); 


/* Adjust for delay introduced by filter */ 

for(i=0; i<num_samples; i++) 
data[i] = data[i+Ml]; 

Psig = calc_power(data, num_samples); 


add_noise(idum, data, num_samples, Psig, EbNo, gain); 

filter(data, M2, Rx, x2, num_samples); 

/* Adjust for delay introduced by filter */ 

for(i=0; i<num_samples; i++) 
data[i] = data[i+M2]; 

demodsoft(data, v, num_bits, soft_metric, Q, a); 

conv_decode(v, conv_n, conv_k, conv_m, conv_num, soft_metric, a); 

deinterleavefv, rows, rs_n, rs_m, I, I_num); 

rs_decode(v, rs_n, rs_k, rs_q, rs_m, rs_num, erasure_counter, rs_counter); 

for (i=0; i<N; i++) 
if (v[i] != u[i]) 

{ 

errors += 1 ; 

} 
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N_total += N; 


outfO=fopen(cFilename,"at+"); 

for (i=0; i<=rs_t; i+=2)fprintf(outf0,"%de= %d ",i, erasure_counter[i]); 

fprintf(outfO,”\nnum = %d\n\nRS counter = ",rs_num); 

for(i=0; i<6; i++) fprintf(outfO," %d",rs_counter[i]); 

fprintf(outfO,"\n\nl. S = 0. No Decoding\n2. deg_lambda > 2*t"); 
fprintf(outf0,"\n3. 2*e + f > 2*t\n4. num_errors != degjambda"); 
fprintf(outf0,"\n5. Correct decoding\n5. Incorrect Decoding"); 


fprintf(outfO,"\n\nEb/No = %2.2f dB",EbNo); 
fprintf(outfO,"\n(%d, %d, %d ) Reed Solomon Code",rs_n, rs_k, rs_t); 
fprintf(outfO,"\ninterleaving depth I = %d",I); 

fprintf(outfO,"\n(%d, %d, %d) convolutional code", conv_n, conv_k, conv_m); 
fprintf(outfO,"\niteration number = %d\n# errors = %d\n N total = %ld", ii+1, errors, N_total); 
Pe = (float )errors/N_total ; 
fprintf(outfO,"\nPe = %2.10An\n",Pe); 
fclose(outfO); 


} 

outfO=fopen(cFilename, " at+" ) ; 
fprintf(outfO,"\nerasure counter = "); 

for (i=0; i<=rs_t; i+=2) fprintf(outfO,"%de= %d ",i, erasure_counter[i]); 

fprintf(outfO,"\nnum = %d\n\nRS counter = ",rs_num); 

for(i=0; i<6; i++) fprintf(outfO," %d",rs_counter[i]); 

fprintf(outfO,"\n\nl. S = 0. No Decoding\n2. deg_lambda > 2*t"); 
fprintf(outf0,"\n3. 2*e + f > 2*t\n4. num_errors !- deg_lambda H ); 
fprintf(outfO, M \n5. Correct deccding\n5. Incorrect Decoding"); 


fprintf(outfO, "\n\nEb/No = %2.2f dB",EbNo); 
fprintf(outfO,"\nNumber of erasures = %d", num_erasures); 
fprintf(outfO,"\n(%d, %d, %d ) Reed Solomon Code",rs_n, rs_k, rs_t); 
fprintf(outfO,"\ninterleaving depth I = %d",I); 

fprintf(outfD,"\n(%d, %d, %d) convolutional code", conv_n, conv_k, conv_m); 
fprintf(outfO,"\niteration number = %d\n# errors = %d\n N total = %ld", ii+1, errors, N_total); 
Pe = (float)errors/N_total; 
fprintf(outfD,"\nPe = %2.10f\n\n",Pe); 

fclose(outfO); 

fclose(inpfl); 

fclose(inpf2); 

fclose(inpf3); 


free(xl); 

free(x2); 

free(data); 
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free(v); 

free(u); 

free(Rx); 

free(Tx); 

free_2d_int_matrix(2, soft_metric); 

} 


B.10 Main program for ideal system using erasure Method 2 

# include <stdio.h> 

# include <math.h> 

# include <stdlib.h> 

# include <stddef.h> 

# include "mmt33.h" 

# define FILE_PRINT 1 /* Prints either to file(l) or screen(O) */ 

/* Note: All DEBUG flags must be 0 */ 

# define GEN.STAT 1 

# define PAUSE 0 /* uses getchar(); to pause the output displayed to the screen */ 

# define L_DEBUG 0 /* prints out reliability values in SOVA */ 

# define DEBUG 0 

# define DEBUG 1 0 

# define DEBUG_DUMP 0 /* Shows the whole deinterleaving frame with reliability values */ 

# define Nss 16 /* Number of samples per symbol */ 

# define Rb 1 /* The bit rate */ 

# define Tb(1.0/Rb) 

# define Eb 1.0 /* The Energy per bit */ 

struct mat 

{ 

int r; 
int r2; 
int reli; 
int flag; 

}; 


/* Defining the random number generator constants */ 

# define IM1 2147483563 

# define IM2 2147483399 

# define AM(1.0/IM1) 

# define IMM1 (IM1-1) 

# define IA1 40014 

# define IA2 40692 

# define IQ 1 53668 

# define IQ2 52774 

# define IR1 12211 

# define IR2 3791 

# define NTAB 32 
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# define NDIV (1+IMM1/NTAB) 

# define EPS 1 .2e-7 

# define RNMX (1.0-EPS) 

int main(int argc, char *argv[]) 

{ 

int *v, *v2, *u, i, ii, j, errors, Q; 

long *idum; 

long idum_cell; 

long num_bits, N, N_total; 

int rs_m, rs_t, rs_n, rs„k, rs_q, rs_num; 

int conv_k, conv_n, convjm, conv_num; 

double Pe, x, var, sd, A, amp, b, a; 

double EbNo, gain, *prob; 

int max_repeat, I, rows, I_num, level, frame_error; 

FILE *outfO; 

int **soft_metric; 

int *erasure_counter, *counter, *dccoded_word; 

int *frame_failure, *error_statl, *error_stat2, **num_bits„per_level; 

char cFilename[80]; 

idum=&idum_cell; 

EbNo=atof(argv[ 1 ]); 

N = atol(argv[2]); 

I = atol(argv[3]); 
max_repeat = atol(argv[4]); 


if(FILE_PRINT) sprintf(cFilename,"B_Method2_%2.3fEbNo%dI.output",EbNo,I); 


/* The Reed Solomon coding parameters *1 

rs_m = 8; 
rs_t = 16; 

rs_q = pow(2, rs_m); 
rs_n = rs_q- 1 ; 
rs_k = rs_n - 2*rs_t; 


/* Interleaving parameters */ 
rows = rs_m*rs_n; 

/* The Convolutional coding parameters */ 

conv_k = 1 ; 
conv_n = 2; 
conv_m = 6; 

Q = 256; /* Number of soft decision levels */ 
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rs_num = N/(rs_k*rs_m); 

if( N % (rs_k*rs_m) != 0) 

{ 

++rs_num; 

N = rs„num*rs_k*rs_m; 

} 


I_num = rs_num/I; 

if(rs_num%I != 0) 

{ 

++ I_num; 
rs_num = I_num*I; 

N = rs_num*rs_k*rs_m; 

} 


gain = (double)(conv_k*rs_k)/(conv„n*rs_n); 
a = 10.0; 

conv_num = (rs_num*rs_n*rs_m +50)/conv_k; 
num_bits = conv_n*(conv_num + 100); 
errors = 0; 

soft_metric = int_matrix_2d(2, Q); 
num_bits_per_level = int_matrix_2d(2, Q); 
counter = ivector(7); 


prob = dvector(Q); 

/* stuff from RS encode *1 
frame_failure = ivector(I+l); 

error_statl = ivector(32); /* Keeps track of how many RSW contain x number of errors */ 
error_stat2 - ivector(32); 

decoded_word = ivector(I); 

for(i=0; i<7; i++) counter[i] = 0; 

amp = 1.0; 

N_total = 0; 

v = ivector(num_bits); 

u = ivector(N); 

v2 = ivector(num_bits); 

*idum = -10; /* Random number initial seed */ 

A = gain*pow(10.0, (EbNo/10.0)); 
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var = 1. 0/(2* A); 


sd = sqrt(var); 

for(ii-0; ii<max_repeat; ii++) 

{ 

for (i=0; i<num_bits; i++) v[i] = 0; 

for (i=0; i<N; i++) 

{ 

v[i] = bitgen(idum); 
u[i] = v[i] ; 

} 

rs_encode(v, rs_n, rs_k, rs_m, rs_num); 

for(i=0; i<num_bits; i++) v2[i] = v[i]; 

interleave^, rows, rs_n, rs_m, I, I_num); 

conv_encode(v, conv_n, conv_k, conv__m, conv__num); 

for (i=0; i<num_bits; i++) 

{ 

if(v[i] == 1) x = 1.0 + sd*gasdev2(idum); 
else x = - 1 .0 + sd*gasdev2(idum); 

/* Demod */ 

for(j=0; j<Q; j++) 

if(x>=amp*(2*j-Q)/Q && x< amp*(2*(j+l)-Q)/Q) 
level =j; 

if(x >= amp) level = Q- 1 ; 
if(x <= -amp) level = 0; 

num_bits_per_Ievel[v[i]][level] += 1; 

v[i] = level; 

} 

for(i=0; i<Q; i++) 

prob[i] = (double)(num_bits_per_level[0][i] + num_bits_per_level[l][Q-i-l] + l)/(2*num_bits +Q) 

b = -log(prob[Q-l])/log(2.0); 

for(i=0; i<Q; i++) 

{ 

soft_metric[0][i] = (int)(floor)(a*(log(prob{i])/log(2.0) + b)); 
soft_metric[l][Q-i-l] = soft_metric[0][i]; 

} 

conv_decode(v, conv_n, conv_k, conv_m, conv_num, soft_metric, a); 
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deinterleave(v, rows, rs_n, rs_m, I, I_num); 

rs_decode_erasure_method_2(v, v2, counter, rs_n, rs_k, rs_q, rs_m, rs_num, 

erasure_counter, I, decoded_word, frame_failure, error_statl, error_stat2); 

for (i=0; i<N; i++) 
if (v[i] != u[i]) 
errors += 1 ; 

N_total += N; 

if(!FILE_PRINT) 

{ 

printf("\n\n"); 

printf(" Erasure Counter\n"); 

for (i=0; i<=2*rs_t-4; i+=2) printf("%de= %d ",i, erasure_counter[i]); 

printf("\n\nError Counter l (number of times x amount of errors occoured)\n"); 

for(i=0; i< 16; i++) 

{ 

if(i==7) printf("\n"); 

printf("%derr = %d ", 16+i, error_statl[i]); 

} 


printf("\n\nError Counter 2 (# of errors for non-decoded word)\n"); 

for(i=0; i< 16; i++) 

{ 

if(i==7) printf("\n"); 

printf("%derr = %d ", 16+i, error_stat2[i]); 

} 


printf("\n\nFrame decoded with x correct\n"); 

for(i=0; i<=I; i++) printf("%dcct = %d ",i, frame„failure[i]); 


if(GEN_STAT) 

{ 

printf("\n\nNumber of GE flagged GE = %d Prob = %1.3f\n", counter[4], 
(float)counter[4]/counter[6] ) ; 

printf("Number of GE flagged BE = %d Prob = % 1 .3f\n", 
counter[5],(float)counter[5]/counter[6]); 

printf("Total number of GE declaired = %d\n",counter[6]); 

printf("\nNumber of BE flagged BE = %d Prob = %1.3f\n", counter[7], 
(float)counter[7]/counter[9]); 

printf("Number of BE flagged GE = %d Prob = % 1 .3f\n", counter[8], 
(float)counter[8]/counter[9]); 

printf("Total number of BE declaired = %d\n",counter[9]); 

} 


printf("\n\nBefore Decoding, no erasures used:\n"); 
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printf("\nTotaI number of bytes = %d", counter[l 1]); 
printf("\nNumber of byte errors = %d", counter[12]); 
printf("\nProb of byte error = %2.51f',(double)counter[12]/counter[l 1]); 

printf("\n\nTotal number of frames = %d", counter[13]); 
printf("\nNumber of frames in error = %d", counter[ 14]); 
printf("\nProb of frame error = %2.51f", (double)counter[14]/counter[13]); 

printf("\n\nNumber of RSW = %d",counter[0]); 

printf("\nNumber of RSW correct = %d",error_statl[0]); 

printf("\nNumber of RSW incorrect = %d", counter[0]-error_statl[0]); 

printf("\nProb of RSW incorrect = %2.5If',(double)(counter[0]-error_statl[0])/counter[0]); 

printf("\n\nAfter decodingAn"); 

frame_error = 0; 

for(i=0; i<I; i++) frame_error += frame_failure[i]; 

printf("\nTotal number of frames = %d", counter[13]); 
printf("\nNumber of frames in error = %d", frame_error); 
printf("\nProb of frame error = %2.51f',(double)frame_error/counter[13]); 


printf("\n\nNumber of RSW - %d",counter[0]); 
printf(”\nNumber of RSW correct = %d",counter[l]); 
printf(’’\nNumber of RSW incorrect = %d", counter[2]); 
printf("\nProb of RSW incorrect = %2.5lf',(double)counter[2]/counter[0]); 

printf("\n\nTotal number of decoding trials = %d",counter[15]); 

printf("\nAverage number of decoding trials = %3.31f", (double)counter[15]/counter[0]); 

printf("\nNumber of decoding failures = %d\n",counter[10]); 


printf("\n\nEb/No = %2.2f dB”,EbNo); 

printf("\n(%d, %d, %d) Reed Solomon Code",rs_n, rs_k, rs_t); 

printf('’\ninterleaving depth I = %d",I); 

printf("\n(%d, %d, %d) convolutional code", conv_n, conv_k, conv_m); 
printf("\niteration number = %d\n# errors = %d\n N total = %ld", ii+1, errors, N_total); 
Pe = (float)errors/N_total; 
printf("\nPe = %2.10f\n\n",Pe); 

} 

if(FILE_PRINT) 

{ 


outfO=fopen(cFilename,"at+"); 
fprintf(outfO, "\n\nErasure C ounter\n " ) ; 

for (i=0; i<=2*rs_t-2; i+=2) fprintf(outfO,"%de= %d ",i, erasure_counter[i]); 
fprintf(outfO,"\n\nError Counter 1 (number of times x amount of errors occoured)\n"); 
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for(i=0; i< 1 6; i++) fprintf(outfO,"%derr = %d ", 16+i, error_statl[i]); 


for(i=0; i< 1 6; i++) 

{ 

if(i==7) fprintf(outfO,"\n"); 
fprintf(outfO,"%derr = %d ", 16+i, error_statl[i]); 

} 


fprintf(outfO,"\n\nError Counter 2 (# of errors for non-decoded word)\n"); 

for(i=0; i< 16; i++) 

{ 

if(i==7) fprintf(outfO,"\n"); 
fprintf(outfO,"%derr = %d ", 16+i, error_stat2[i]); 

} 


fprintf(outfO,"\n\nFrame decoded with x correct\n"); 

for(i=0; i<=I; i++) fprintf(outfO,"%dcct = %d ",i, frame_failure[i]); 

if(GEN_STAT) 

{ 

fprintf(outfO,"\n\nNumber of GE flagged GE = %d Prob = % 1 .3f\n", counter[4], 
(float)counter[4]/counter[6]); 

fprintf(outfO, "Number of GE flagged BE = %d Prob = %1.3f\n", 
counter[5],(float)counter[5]/counter[6]); 

fprintf(outfO, "Total number of GE declaired = %d\n",counter[6]); 

fprintf(outfO,"\nNumber of BE flagged BE = %d Prob = %1.3f\n", counter[7], 
(float)counter[7]/counter[9j); 

fprintf(outfO, "Number of BE flagged GE = %d Prob = % 1 ,3f\n", counter[8], 
(float)counter[8]/counter[9]); 

fprintf(outfO, "Total number of GE declaired = %d\n",counter[9]); 

} • • 


fprintf(outfO,"\n\nBefore Decoding, no erasures used:\n"); 

fprintf(outfO,"\nTotal number of bytes = %d", counter[l 1]); 
fprintf(outfO,"\nNumber of byte errors = %d", counter[12]); 
fprintf(outfO,"\nProb of byte error = %2.51f',(double)counter[12]/counter[l l]); 

fprintf(outfO,"\n\nTotal number of frames = %d", counter[13j); 
fprintf(outfO,"\nNumber of frames in error = %d", counter[14]); 
fprintf(outfO,"\nProb of frame error = %2.51f, (double)counter[14]/counter[13]); 

fprintf(outfO,"\n\nNumber of RSW = %d",counter[0]); 
fprintf(outfO,"\nNumber of RSW correct = %d",error_statl [0]); 
fprintf(outfO,"\nNumber of RSW incorrect = %d", counter[0]-error_statl[0]); 
fprintf(outfO,"\nProb of RSW incorrect = %2.51f',(double)(counter[0]-error_statl[0])/counter[0]); 

fprintf(outfO, "\n\n After decoding :\n " ) ; 

frame_error - 0; 
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for(i=0; i<I; i++) frame_error += frame_failure[i]; 

fprintf(outfO,"\nTotal number of frames = %d", counter[13]); 
fprintf(outfO,"\nNumber of frames in error = %d", frame_error); 
fprintf(outfO,"\nProb of frame error = %2.51f’,(double)frame_error/counter[13]); 


fprintf(outfO,"\n\nNumber of RSW = %d",counter[0]); 
fprintf(outfO,"\nNumber of RSW correct - %d",counter[l]); 
fprintf(outfO,"\nNumber of RSW incorrect = %d", counter[2j); 
fprintf(outfO,"\nProb of RSW incorrect = %2.51f",(double)counter[2J/counter[0]); 

fprintf(outfO,"\n\nTotal number of decoding trials = %d",counter[15]); 
fprintf(outfO,"\nAverage number of decoding trials = %3.31f, (double)counter[15]/counter[0]); 

fprintf(outfO,"\nNumber of decoding failures = %d\n", counter! 10]); 


fprintf(outfO,"\n\nEb/No = %2.2f dB",EbNo); 
fprintf(outfO,"\n(%d, %d, %d) Reed Solomon Code",rs_n, rs_k, rs_t); 
fprintf(outfO,"\ninterleaving depth I = %d",I); 

fprintf(outfO,"\n(%d, %d, %d) convolutional code", conv_n, conv_k, conv_m); 
fprintf(outfO,"\niteration number = %d\n# errors = %d\n N total = %ld", ii+1, errors, N_total); 
Pe = (float)errors/N_total; 
fprintf(outfO,"\nPe = %2. 10f\n\n",Pe); 

fclose(outfO); 

} 


free(v); 

free(u); 

free(v2); 

free_2d_int_matrix(2, soft_metric); 

free_2d_int_matrix(2, num_bits_per_level); 

free(decoded_word) ; 

free(frame_failure); 

free(error_statl); 

free(error_stat2); 

} 
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